# Un cadre général pour le traitement des problèmes d'interface

Hachmi Ben Dhia, Marc Kham, Patrick Massin, Mohamed Torkhani

► To cite this version:

Hachmi Ben Dhia, Marc Kham, Patrick Massin, Mohamed Torkhani. Un cadre général pour le traitement des problèmes d'interface. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421846

HAL Id: hal-01421846

https://hal.science/hal-01421846v1

Submitted on 23 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Un cadre général pour le traitement des problèmes d'interface

H. Ben Dhia¹, M. Kham², P. Massin³, M. Torkhani¹,³

¹ Ecole Centrale Paris, LMSSMAT

Grande Voie des Vignes, 92295 Châtenay-Malabry Cedex
bendhia@mssmat.ecp.fr

² Département Analyses Mécaniques et Acoustiques, EDF R&D
1, avenue du général de Gaulle, 92141 Clamart Cedex
marc.kham@edf.fr

³ Laboratoire de Mécanique des Structures Industrielles Durables, UMR EDF/CNRS 2832
1, avenue du général de Gaulle, 92141 Clamart Cedex
{ patrick.massin, mohamed.torkhani } @edf.fr

**Résumé** – L’enjeu de ce travail est de proposer un cadre général pour le traitement des problèmes d’interface. Le problème de contact-frottement est écrit de manière continue, puis exprimé sous forme variationnelle pour un traitement de type élément fini. Les non linéarités de contact-frottement dues à la géométrie, au statut des points de contact et au seuil de frottement sont découplées. L’approche est utilisée à la fois pour le traitement de problèmes de contact avec ou sans usure, avec ou sans dynamique, en FEM ou en X-FEM.

**Mots clés** – X-FEM, contact frottant, grands glissements, usure, 3D, condition LBB.

## 1 Introduction

Le traitement du contact reste toujours un défi pour grand nombre de codes numériques, soit en terme de qualité des résultats, soit en terme de capacité de calculs. Dans le contexte d'une collaboration entre ECP et EDF R&D qui dure depuis maintenant une dizaine d'années, une discrétisation d'une famille de formulations lagrangiennes continues du problème mécanique de contact frottant, englobant la formulation lagrangienne augmentée [3], a été implantée dans Code_Aster, Code de Mécanique de la R&D d'EDF.

Cette famille lagrangienne est dérivée d'une formulation variationnelle des lois fortes de contact frottant, écrites sous forme d'équations locales, grâce à l'introduction de champs de signes (ou champs de niveaux) qualifiant les différents statuts (contactant/non contactant et glissant/adhérant) de chaque point des surfaces potentielles de contact [1, 9].

Les problèmes discrets sont obtenus par usage de différentes méthodes d'approximation. Ainsi, des éléments finis compatibles sont utilisés pour l'approximation des champs de déplacement (ou de vitesse) et de densités d'efforts de
contact. La méthode des différences finies est utilisée pour les discrétisations temporelles (vitesse) et une méthode de collocation, cohérente avec une méthode d'intégration numérique précise des termes énergétiques de contact, est retenue pour l'approximation des champs (irréguliers) de signes [1, 9].

Les principales étapes de l'algorithme numérique utilisé pour la résolution des problèmes discrets non linéaires et non réguliers, ainsi obtenus, peuvent être présentées comme quatre boucles imbriquées à chaque pas de temps :

* boucle de réactualisation de la géométrie des surfaces de contact et lancement de l'algorithme d'appariement;

* boucle sur les seuils de frottement (méthode de point-fixe);

* boucle sur les statuts de contact (méthode des contraintes actives [2]);

* boucle de Newton généralisée.

En scindant les non linéarités de contact, cette approche vise la robustesse. Par ailleurs, elle génère des matrices tangentes symétriques.

Des travaux récents [4] ont permis de prendre en compte l'usure associée au frottement glissant, en introduisant une extension vectorielle locale de la loi phénoménologique scalaire d'Archard [5] au cas des grandes transformations, des grands glissements et des grandes usures. En statique, la formulation a pu être étendue aux éléments finis X-FEM, avec petits déplacements [6] ou grandes transformations [7], même si des différences essentielles apparaissent dans le choix d'approximation des espaces d'éléments finis pour satisfaire la condition LBB [8].

En dynamique, deux formulations des lois de contact peuvent être proposées : l'une basée sur les déplacements (Signorini) et l'autre sur les vitesses (Moreau) (cf. [9]), plus adaptée que la précédente à la dynamique du contact (chocs entre solides déformables).

## 2 Problème de contact frottant et formulation variationnelle mixte associée

Nous nous plaçons dans le cadre des petites déformations, petites usures, et les solides en contact sont considérés, dans les applications présentées ici, comme élastiques linéaires. Une loi de coulomb est utilisée pour modéliser le frottement. Les déplacements **u** peuvent être grands, notamment au niveau de l'interface de contact. Une approche de type maître-esclave est utilisée pour le traitement des surfaces de contact. On tient compte du champ d'usure normal à la surface de contact $w_n$. Nous notons $\mathbf{x}_w^1$ et $\bar{\mathbf{x}}_w^1$ les positions courantes usées d'un point $\mathbf{p}^1$ de la surface de référence du solide esclave et de son vis-à-vis apparié $\bar{\mathbf{p}}^1$ de la surface de référence du solide maître, telles que :

$$ \mathbf{x}_w^1 = \mathbf{x}^1(\mathbf{p}^1, t) - w_n^1(\mathbf{p}^1, t)\mathbf{n}(\mathbf{p}^1, t) \text{ et } \bar{\mathbf{x}}_w^1 = \bar{\mathbf{x}}^1(\bar{\mathbf{p}}^1, t) - w_n^2(\bar{\mathbf{p}}^1, t)\mathbf{n}(\bar{\mathbf{p}}^1, t) \text{ où } \mathbf{n}(\mathbf{p}^1, t) \text{ et } \mathbf{n}(\bar{\mathbf{p}}^1, t) \text{ sont les normales sortantes aux solides aux points } \mathbf{p}^1 \text{ et } \bar{\mathbf{p}}^1. $$

On définit alors la distance $d^w$ entre les points appariés $\mathbf{p}^1$ et $\bar{\mathbf{p}}^1$ sur les corps usés par $d^w = (\mathbf{x}_w^1 - \bar{\mathbf{x}}_w^1) = [[\mathbf{x}]] = [(\mathbf{x}_w^1 - \bar{\mathbf{x}}_w^1)\mathbf{n}]\mathbf{n} + \mathbf{d}_v$, avec $\mathbf{n} = -\mathbf{n}(\bar{\mathbf{p}}^1, t)$ normale unitaire au solide maître au point $\bar{\mathbf{p}}^1$ dirigée vers l'intérieur du solide maître. La distance $d^w$ se réduit à
$\mathbf{d}^w = (\mathbf{x}_w^1 - \bar{\mathbf{x}}_w^1) = (d_n - (w_n^1 + w_n^2))\mathbf{n}$ en cas de contact effectif, où l'on note
$d_n^w = (d_n - (w_n^1 + w_n^2)) = d_n - w_n$ le jeu normal usé. Un modèle d'usure type Archard peut être
utilisé ici, mais des modèles plus évolués en transformations finies sont présentés dans
[4].

Dans le cadre élément fini qui nous intéresse, les équations du contact au niveau des surfaces
potentielles de contact sont écrites sous forme variationnelle, en faisant intervenir un
multiplicateur de Lagrange $\lambda$ pour la pression de contact et un vecteur semi-multiplicateur de
frottement $\Lambda$ de norme comprise entre 0 et 1 (0 en non contactant, 1 en glissement, entre 0 et 1
en adhérence) et dont la direction est celle de la réaction tangentielle de frottement (la valeur de
la pression tangentielle est donnée par le produit de la pression de contact par la norme du semi-
multiplicateur). En quasi-statique, le problème à résoudre se ramène finalement à trouver les
champs $(\mathbf{u}_k, \lambda_k, \mathbf{\Lambda}_k, w_{nk}) \in \text{CA} \times H_c \times \mathbf{H} \times H_w$ tels que
$$ \forall (\mathbf{u}_k^*, \lambda_k^*, \mathbf{\Lambda}_k^*, w_{nk}^*) \in \text{CA} \times H_c \times \mathbf{H} \times H_w : $$

$$
\begin{align}
\mathbf{G}_{\text{int}}(\mathbf{u}_k, \mathbf{u}^*) + \mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) + \mathbf{G}_{\text{frot}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) &= 0 \tag{1} \\
\mathbf{G}_{\text{cont}}^{\text{faible}}(\lambda_k, \mathbf{u}_k, w_{nk}, \lambda^*) &= 0 \nonumber \\
\mathbf{G}_{\text{frot}}^{\text{faible}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{\Lambda}^*) &= 0 \nonumber
\end{align}
$$

où l'on reconnaît le système d'équilibre, prenant en compte les contributions de
contact frottement et l'usure, et la forme faible des équations de contact et de frottement.

On donne ci-dessous les expressions de ces différents termes :

$$
\begin{align}
\mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) &= -\int_{\Gamma_c} S_{uk}^w \cdot g_{nk}^w \cdot \llbracket \mathbf{u}^* \rrbracket_{\nu} \cdot d\Gamma_c \tag{2} \\
\mathbf{G}_{\text{frot}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) &= -\int_{\Gamma_c} S_{uk}^w \cdot \mu \cdot \lambda_k \cdot \left( S_{fk} \cdot \mathbf{\Lambda}_k + (1 - S_{fk}) \cdot \frac{\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{fk}}{\| \mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{fk} \|} \right) \cdot \llbracket \mathbf{u}^* \rrbracket_{\tau} \cdot d\Gamma_c \nonumber \\
\mathbf{G}_{\text{cont}}^{\text{faible}}(\lambda_k, \mathbf{u}_k, w_{nk}, \lambda^*) &= -\frac{1}{h_n} \int_{\Gamma_c} (\lambda_k - S_{uk}^w \cdot (\lambda_k - h_n \cdot d_{nk}^w)) \lambda^* \cdot d\Gamma_c \nonumber \\
\mathbf{G}_{\text{frot}}^{\text{faible}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{\Lambda}^*) &= \frac{1}{h_\tau} \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk}^w \cdot \mathbf{\Lambda}_k \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c - \frac{1}{h_\tau} \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk}^w \cdot S_{fk} \cdot (\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{fk}) \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c \nonumber \\
&\quad - \frac{1}{h_\tau} \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk}^w \cdot (1 - S_{fk}) \cdot \frac{\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{fk}}{\| \mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{fk} \|} \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c + \int_{\Gamma_c} (1 - S_{uk}^w) \cdot \mathbf{\Lambda}_k \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c \nonumber
\end{align}
$$

où $h_n$ est un réel quelconque, $g_{nk}^w = \lambda_k - \rho_n \cdot d_{nk}^w$ est le multiplicateur de contact augmenté et
$\mathbf{v}_{fk} = \Delta[\llbracket \mathbf{x}_k \rrbracket_{\tau} / \Delta t_k$ est la vitesse relative de glissement entre les solides maître et esclave, avec
$\mathbf{x}_{\tau} = (\mathbf{1} - \mathbf{n} \otimes \mathbf{n}) \cdot \mathbf{x}$. On introduit en outre les deux fonctions indicatrices suivantes :

$$
S_{uk}^{w} = \chi(g_{nk}^{w}) = \begin{cases} 1 & g_{nk}^{w} \le 0 \text{ contact} \\ 0 & g_{nk}^{w} > 0 \text{ pas de contact} \end{cases} \tag{3}
$$

avec $\rho_n > 0$ paramètre réel d'homogénéisation et :

$$
S_{fk} = \chi(\|\mathbf{\Lambda}_k + \rho_\tau \cdot \mathbf{v}_{fk}\| - 1) = \begin{cases} 1 & \|\mathbf{\Lambda}_k + \rho_\tau \cdot \mathbf{v}_{fk}\| \le 1 \text{ adhérent} \\ 0 & \|\mathbf{\Lambda}_k + \rho_\tau \cdot \mathbf{v}_{fk}\| > 1 \text{ glissant} \end{cases} \quad (4)
$$
avec $\rho_{\tau} > 0$ paramètre réel d'homogénéisation.

Cette formulation peut être modifiée pour prendre en compte une interface représentée par une couche d'aspérités d'épaisseur a. Les deux solides maître et esclave entrent alors en contact dès que le jeu devient inférieur à l'épaisseur des aspérités. Tant que le jeu n'est pas nul, l'écrasement de la couche d'aspérités est modélisé par un modèle de compliance. Le terme $\mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*)$ est alors modifié de la manière suivante :

$$ \mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) = -\int_{\Gamma_c} S_{uk}^{w} g_{nk}^{w} [[\mathbf{u}^*]]_n \cdot d\Gamma_c + \int_{\Gamma_c} (S_{uk} \cdot \kappa_n \cdot (d_{nk} + a)^2) [[\mathbf{u}^*]]_n \cdot d\Gamma_c + \int_{\Gamma_c} (S_{uk} \cdot c_n \cdot (d_{nk} + a)) [[\mathbf{v}]]_n \cdot d\Gamma_c + \int_{\Gamma_c} (S_{uk} \cdot c_n \cdot (d_{nk} + a)) [[\mathbf{u}^*]]_n \cdot d\Gamma_c \quad (5) $$

avec $S_{uk} = \chi(d_{nk} + a)$. Les paramètres $(\kappa_n, c_n)$, strictement positifs, peuvent dans certains cas être identifiés expérimentalement, le paramètre a étant lié à la hauteur moyenne des aspérités.

La formulation peut être aussi stabilisée par le rajout d'un terme de type pénalisation (dont le paramètre est déconnecté de celui d'augmentation) au niveau de $\mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*)$ et $\mathbf{G}_{\text{frot}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*)$ apparaissant dans les équations d'équilibre. Elle donne alors lieu à une formulation lagrangienne stabilisée [9]; le terme de stabilisation n'introduisant aucune perturbation de la solution. Si on ne garde que les termes de stabilisation, on retrouve une formulation par pénalisation classique, qui n'intervient qu'au niveau de l'équation d'équilibre : dans ce cas là on retrouve la dépendance des résultats obtenus au choix des coefficients de pénalisation $\kappa_n$ et $\kappa_{\tau}$ strictement positifs. Les termes $\mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*)$ et $\mathbf{G}_{\text{frot}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*)$ sont alors modifiés de la manière suivante :

$$ \begin{gathered} \mathbf{G}_{\text{cont}}(\lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) = -\int_{\Gamma_c} S_{uk}^{w} g_{nk}^{w} [[\mathbf{u}^*]]_n \cdot d\Gamma_c + \int_{\Gamma_c} (S_{uk}^{w} \cdot \kappa_n \cdot d_{nk}) [[\mathbf{u}^*]]_n \cdot d\Gamma_c \quad (6) \\ \mathbf{G}_{\text{frot}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{u}_k, w_{nk}, \mathbf{u}^*) = -\int_{\Gamma_c} S_{uk}^{w} \cdot \mu \cdot \lambda_k \cdot \left( S_{jk} \cdot \mathbf{\Lambda}_k + (1 - S_{jk}) \frac{\mathbf{\Lambda}_k + h_{\tau} \cdot \mathbf{v}_{jk}}{\|\mathbf{\Lambda}_k + h_{\tau} \cdot \mathbf{v}_{jk}\|} \right) [[\mathbf{u}^*]]_{\tau} \cdot d\Gamma_c - \int_{\Gamma_c} S_{uk}^{w} S_{jk} \cdot \kappa_{\tau} \cdot \mu \cdot \lambda_k \cdot \mathbf{v}_{jk} [[\mathbf{u}^*]]_{\tau} \cdot d\Gamma_c \end{gathered} $$

En dynamique, en l'absence d'usure, on modifie le formalisme précédent pour prendre en compte les lois de Signorini-Moreau introduisant des conditions sur la vitesse de contact. La forme variationnelle des équations d'équilibre et de contact-frottement devient alors :

$$ \begin{align} \mathbf{G}_{\text{cont}}^{\text{dyna}}(\lambda_k, \mathbf{v}_k, \mathbf{v}^*) &= -\int_{\Gamma_c} S_{uk} \cdot S_{vk} \cdot \lambda_k \cdot [[\mathbf{v}^*]]_n \cdot d\Gamma_c \tag{7} \\ \mathbf{G}_{\text{frot}}^{\text{dyna}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{v}_k, \mathbf{v}^*) &= \int_{\Gamma_c} S_{uk} \cdot \mu \cdot \lambda_k \cdot \left( S_{jk} \cdot \mathbf{\Lambda}_k + (1 - S_{jk}) \frac{\mathbf{g}_{jk}}{\|\mathbf{g}_{jk}\|} \right) [[\mathbf{v}^*]]_\tau \cdot d\Gamma_c \nonumber \\ \mathbf{G}_{\text{cont}}^{\text{dyna,faible}}(\lambda_k, \mathbf{v}_k, \lambda^*) &= -\frac{1}{h_\tau} \int_{\Gamma_c} (\lambda_k - S_{uk} \cdot S_{vk} \cdot (\lambda_k - h_n \cdot [[\mathbf{v}_k]]_n)) \lambda^* \cdot d\Gamma_c \nonumber \\ \mathbf{G}_{\text{frot}}^{\text{dyna,faible}}(\mathbf{\Lambda}_k, \lambda_k, \mathbf{v}_k, \mathbf{\Lambda}^*) &= \frac{1}{h_\tau} \cdot \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk} \cdot \mathbf{\Lambda}_k \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c - \frac{1}{h_\tau} \cdot \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk} \cdot S_{jk} \cdot (\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{jk}) \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c \nonumber \\ &\quad - \frac{1}{h_\tau} \cdot \int_{\Gamma_c} \mu \cdot \lambda_k \cdot S_{uk} \cdot (1 - S_{jk}) \cdot \frac{\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{jk}}{\|\mathbf{\Lambda}_k + h_\tau \cdot \mathbf{v}_{jk}\|} \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c + \int_{\Gamma_c} (1 - S_{uk}) \cdot \mathbf{\Lambda}_k \cdot \mathbf{\Lambda}^* \cdot d\Gamma_c \nonumber \end{align} $$

où $S_{uk} = \chi(-d_{nk})$ et $S_{vk} = \chi(\lambda - \rho_n \cdot [[\mathbf{v}_k]]_n)$. Des résultats sont donnés dans une communication du présent Colloque intitulée « Formulation continue en vitesse du contact-frottement ».

L'adaptation de cette formulation continue [1] au cadre X-FEM est relativement aisée et développée dans la communication intitulée « Modélisation de discontinuités 3D avec contact-
frottement en grandes transformations dans le cadre X-FEM ». On ne la rappelle donc pas ici et
on renvoie aux articles [6] et [7].

3 Résolution numérique

La résolution des problèmes formulés ci-dessus nécessite des méthodes de discrétisation en espace et éventuellement en temps, ainsi que des algorithmes itératifs.

Pour le premier aspect, nous renvoyons à [9] et aux références qui sont mentionnées dans ce papier. Soulignons toutefois un aspect essentiel soulevé et traité dans [1], lié à la nécessité d'intégration précise des termes de contact pour le passage du pach test de Taylor. Soulignons également qu'en dynamique, divers schémas numériques en temps ont été utilisés et que les meilleurs résultats ont été obtenus avec des schémas d'ordre 1 (θ schéma) ou un schéma de la famille de Newmark intégrant les travaux réalisés dans [10].

Concernant l’aspect algorithmique, la méconnaissance des vis-à-vis en contact est traitée par le biais d’une boucle de point fixe sur la géométrie, combinée à l’usage d’algorithme locaux de Newton. Le caractère non différentiable lié au contact unilatéral est résolu via une méthode de type contraintes actives [2]. Les non linéarités de frottement sont traitées par un point fixe sur le seuil de frottement, ramenant à des modèles de frottement de type Tresca. Le caractère non différentiable lié à la projection sur la boule unité du semi_multiplicateur de frottement est traité par une méthode de type Newton généralisé ou module tangent.

4 Conclusions

Un cadre général et évolutif a été présenté pour le traitement des problèmes de contact-
frottement d’interfaces. De nombreux tests de validation élémentaires sont consultables sur
[www code-aster org](http://www.code-aster.org) ou publiés dans [1][4][6][7] et [9]. Ils montrent que l’on a pu valider la
formulation à la fois en statique et en dynamique, l’étendre aux éléments finis de structures, aux
éléments X-FEM et au traitement de l’usure avec les adaptations nécessaires, preuve de la
souplesse de cette formulation introduite dans un environnement de code industriel.

Références

[1] H. Ben Dhia and M. Zarroug, “Hybrid frictional contact particle-in elements”, *Revue Européenne des Éléments Finis*, vol. 9, 417-430, 2002.

[2] G. Dumont, “ Algorithme des contraintes actives et contact unilatéral sans frottement”,
*Revue Européenne des Éléments Finis*, vol. 4, n°1, 55-73, 2001.

[3] P. Alart and A. Curnier, “A mixed formulation for frictional contact problems prone to Newton like solution methods”, *Comp. Meth. Appl. Meth. Engng.*, Vol. 92, 353-375, 1991.

[4] H. Ben Dhia and M. Torkhani, “Multiscale-Multimodel analyses and wear of thin structures in contact problems”, submitted.

[5] J.F. Archard, "Contact and rubbing of flat surfaces", *J. Appl. Phys.*, Vol. 24, (1953) 981-988.

[6] S. Geniaut, P. Massin and N. Moës, “A stable 3D contact formulation for cracks using X-FEM”, *Revue Européenne de Mécanique Numérique*, vol. 16, n°2, 259-275, 2007.

[7] Nistor I., Guiton M.L.E., Massin P., Moës N. and Géniaut S., “An X-FEM approach for large sliding contact along discontinuities”, accepted in *Int. J. Numer. Meth. Engng.*, 2008.
[8] J.M. Solberg, P. Papadopoulos : "An analysis of dual formulations for the finite element solution of two-body contact problems", *Computational Methods and Applied Mechanics Engineering*, vol. 194, pp. 2734-2780, 2005.

[9] H. Ben Dhia, C. Zammali, "Level-Sets fields, placement and velocity based formulations of contact-impact problems", *International Journal for Numerical Methods in Engineering*, Vol. 69, 2711-2735, 2007.

[10] C. Hager, S. Hüeber, B. Wohlmuth, "A stable energy conserving approach for frictional contact problems based on quadrature formulas", *Internat. J. Numer. Methods Engrg.*, 73, 205-225, 2008.