# Un algorithme efficace pour les problèmes d'impact avec frottement

Zhi-Qiang Feng, Jean-Michel Cros, Benoit Magnain

► To cite this version:

Zhi-Qiang Feng, Jean-Michel Cros, Benoit Magnain. Un algorithme efficace pour les problèmes d'impact avec frottement. Revue Européenne des Éléments Finis, 2005, 14 (1), pp.65-86. 10.3166/reef.14.65-86. hal-01179492

HAL Id: hal-01179492

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

Submitted on 28 Jun 2018

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 algorithme efficace pour les problèmes d'impact avec frottement

Zhi-Qiang Feng — Jean-Michel Cros — Benoît Magnain

Laboratoire de Mécanique et d'Énergétique d'Évry

Université d'Évry-Val d'Essonne

40 rue du Pelvoux

F-91020 Évry

{feng,cros,magnain}@iup.univ-evry.fr

RÉSUMÉ. La méthode du bipotentiel a été employée avec succès pour modéliser des problèmes de contact avec frottement en statique. Cet article présente l'extension de cette méthode pour l'analyse dynamique des problèmes d'impact en présence de frottement. Un algorithme du premier ordre est utilisé pour l'intégration temporelle de l'équation discrétisée du mouvement. Les résultats numériques montrent que l'algorithme préserve quasi parfaitement le principe de conversation de l'énergie et cela sans utiliser aucune technique de régularisation. De plus, il est possible de quantifier la dissipation physique d'énergie sous l'effet du frottement entre les solides en contact.

ABSTRACT. The bipotential method has been successfully applied for the modelling of frictional contact problems in static cases. This paper presents the application of this method for dynamic analysis of impact problems with friction. Instead of second order algorithms, a first order algorithm is applied for the numerical integration of the time-discretized equation of motion. The numerical results prove that the algorithm preserves quasi perfectly the principle of energy conservation without any regularization. In addition, it is possible to quantify the physical energy dissipation introduced by frictional effects between the solids in contact.

MOTS-CLÉS : impact, dissipation d'énergie, méthode du bipotentiel, intégration temporelle.

KEYWORDS: impact, energy dissipation, bipotential method, time-integration.
# 1. Introduction

L'analyse des problèmes de contact avec frottement est très importante pour beaucoup d'applications industrielles. Au cours des deux dernières décennies, des avancées importantes ont été accomplies dans l'analyse des problèmes de contact par la méthode des éléments finis. Un grand nombre d'algorithmes de résolution a été présenté dans la littérature. Le lecteur trouvera des exposés de synthèse parmi les références bibliographiques [KLA 93, WRI 95], ainsi que des monographies par Kikuchi et Oden [KIK 88], Zhong [ZHO 93] et Wriggers [WRI 02]. Les méthodes de pénalisation ou de régularisation [CHA 86, PAR 89] semblent, à première vue, les plus appropriées pour traiter les problèmes de contact. Mais dans ce type de méthode, les conditions de contact et les lois de frottement ne sont pas satisfaites exactement. De plus, il est délicat pour les utilisateurs de choisir le facteur de pénalité adéquat. Les expériences ont montré que ces méthodes présentent des inconvénients en ce qui concerne la stabilité et la précision numérique, en particulier pour tout ce qui touche à la simulation des phénomènes de frottement. Pour pallier ces insuffisances, une méthode du Lagrangien augmenté a été développée par Alart et Curnier [ALA 91]. Cette méthode consiste à déterminer les inconnues (déplacements et réactions) simultanément en utilisant un algorithme de Newton généralisé. Simo et Laursen [SIM 92] ont également proposé une méthode similaire. De Saxcé et Feng ont proposé une méthode du bipotentiel fondée sur la théorie MSI (Matériau Standard Implicite) [DES 91], dans laquelle une nouvelle formulation du Lagrangien augmenté est développée. Pour les problèmes de contact unilatéral avec frottement, la méthode du bipotentiel n'utilise qu'un seul principe variationnel sur le déplacement et une seule inégalité sur le contact unilatéral et le frottement. Dans la méthode du bipotentiel, le problème de contact avec frottement est traité dans un système réduit par un algorithme d'Uzawa avec une seule phase de prédiction-correction sur le cône de frottement [DES 98]. Pour plus de détails sur la formulation du Lagrangien augmenté, nous renvoyons vers les articles de Klarbring et al. [KLA 92, CHR 98].

Pour la résolution des problèmes de dynamique des structures, les schémas d'intégration du second ordre sont les plus utilisés. Parmi ces schémas, on peut citer les algorithmes de Newmark, Wilson et HHT. Hughes et al. [HUG 76] ont proposé une version modifiée de l'algorithme de Newmark pour traiter les problèmes de contact et d'impact. Wriggers et al. [WRI 90] ont développé un algorithme de retour-radial pour traiter les problèmes d'impact. Armero et Petocz [ARM 98] et Laursen et Chawla [LAU 97] ont traité des problèmes dynamiques d'impact sans frottement en tenant compte de la conservation de l'énergie. Plus récemment, Laursen et Love [LAU 02] ont proposé un schéma d'intégration implicite, dans lequel la vitesse est corrigée de sorte à éviter l'interpénetration des surfaces de contact. Dans tous les travaux cités précédemment, les auteurs emploient des schémas d'intégration du second ordre dans lesquels on suppose que l'accélération est constante ou varie de façon linéaire. Or, pour les problèmes d'impact, on constate que la vitesse et l'accélération ne sont pas continues en raison du changement brutal des conditions de contact (impact, rebond,...). Dans ces conditions, les algorithmes du second ordre peuvent conduire à de graves
erreurs. Afin de contourner ces défauts, plusieurs algorithmes du premier ordre ont été proposés par Zienkiewicz et al. [ZIE 84] et Jean [JEA 99].

L'objectif de cet article est de traiter les problèmes d'impact dans le cadre de la dynamique non régulière. À cette fin, on propose un algorithme qui combine la méthode du bipotentiel et un schéma d'intégration temporelle du premier ordre. L'algorithme a été implanté dans le code de calcul par élément finis FER/Impact, ce dernier est écrit en C++ et utilise pleinement les techniques de programmation orientée objets. Deux exemples sont étudiés afin de valider les méthodes proposées. Le premier exemple traite de l'impact oblique avec rebond d'une plaque élastique sur une surface rigide. Le second exemple simule l'impact d'un disque entre deux plans rigides. Afin de rendre compte de la dissipation d'énergie par frottement, les deux exemples sont traités en considérant les contacts avec et sans frottement.

## 2. Description du problème

### 2.1. Cinématique du contact

Dans la suite, les notations et définitions sont introduites. On considère deux solides $B^\alpha$ (figure 1), avec $\alpha = 1, 2$. Chacun d'eux occupe un domaine borné $\Omega^\alpha \subset \mathbb{R}^3$, et un point du domaine est noté $X^\alpha$. En outre, ces solides sont considérés comme élastiques et sont soumis à des grands déplacements. La frontière $\Gamma^\alpha$ de chaque corps est supposée suffisamment régulière de sorte à pouvoir définir un vecteur normal unitaire, noté $n^\alpha$, en chaque point $M$ de $\Gamma^\alpha$. On note $\mathbf{I} = [0, T]$ l'intervalle de temps au cours duquel le chargement est appliqué. À chaque instant, $t \in \mathbf{I}$, la frontière $\Gamma^\alpha$ du solide $B^\alpha$ peut, en général, être découpée en trois parties distinctes : $\Gamma_u^\alpha, \Gamma_t^\alpha, \Gamma_c^\alpha$. Avec, $\Gamma_u^\alpha$ la partie où est imposé le déplacement $\bar{u}^\alpha$, $\Gamma_t^\alpha$ la partie où est imposé l'effort $\bar{t}^\alpha$, et enfin les surfaces de contact $\Gamma_c^\alpha$ où les deux corps $B^1$ et $B^2$ peuvent être en contact à l'instant $t$ :

$$ \Gamma^{\alpha} = \Gamma_{u}^{\alpha} \cup \Gamma_{t}^{\alpha} \cup \Gamma_{c}^{\alpha} \quad [1] $$

Les configurations déformées successives de $B^\alpha$ sont décrites à chaque instant par le champ de déplacement $\mathbf{u}^\alpha$ défini sur $\bar{\Omega}^\alpha$ (c'est-à-dire la fermeture de $\Omega^\alpha$). Sur la surface de contact, une normale unique, notée $\mathbf{n}$, dirigée vers $B^1$ ($\mathbf{n} \equiv \mathbf{n}^2$) est définie et le plan tangent, orthogonal à $\mathbf{n}$ dans $\mathbb{R}^3$, est noté $\mathbf{T}$. Pour construire une base locale orthonormée, deux vecteurs notés $\mathbf{t}_x$ and $\mathbf{t}_y$ sont définis dans le plan $\mathbf{T}$. Pour décrire le contact avec frottement qui peut se produire sur $\Gamma_c$, nous introduisons la vitesse relative par rapport à $B^2$

$$ \dot{\mathbf{u}} = \dot{\mathbf{u}}^1 - \dot{\mathbf{u}}^2 \quad [2] $$

où $\dot{\mathbf{u}}^1$ et $\dot{\mathbf{u}}^2$ sont respectivement les vitesses instantanées de $B^1$ et $B^2$. Soit $\mathbf{r}$ la distribution des forces de contact exercée sur $B^1$ à $M$ par $B^2$. En accord avec le principe d'action-réaction, $B^2$ est soumis au vecteur contrainte $-\mathbf{r}$. Dans le système de coordonnées
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Cinématique du contact</figcaption></figure>

données local défini par le plan tangent **T** et la normale **n**, les vecteurs **u** et **r** peuvent
alors être décomposés de façon unique comme suit :

$$
\begin{align*}
\dot{\mathbf{u}} &= \dot{\mathbf{u}}_t + \dot{u}_n \mathbf{n}, & \dot{\mathbf{u}}_t \in \mathbf{T}, & \dot{u}_n \in \mathbb{R} \\
\mathbf{r} &= \mathbf{r}_t + r_n \mathbf{n}, & \mathbf{r}_t \in \mathbf{T}, & r_n \in \mathbb{R}
\end{align*}
\tag{3}
$$

**2.2. Loi de contact et modèle de frottement**

La loi de contact unilatéral est caractérisée par une condition géométrique de non-
pénétration, un état statique de non-adhésion et un état mécanique complémentaire.
Ces trois conditions sont connues sous le nom des conditions de Signorini. La condi-
tion de non-pénétration contraint le champ de déplacement, elle est donnée par :

$$
g(\mathbf{X}^1) = (\mathbf{X}^1 - \mathbf{X}^2) \cdot \mathbf{n} \geq 0 \quad [4]
$$

où

$$
\mathbf{X}^{\alpha}(t) = \mathbf{X}^{\alpha}(t=0) + \mathbf{u}^{\alpha} \quad [5]
$$

Le point **X**² est défini comme la projection du point **X**¹ ∈ Γ¹_c sur la surface Γ²_c. En
notant h l'écart initial :

$$
h = (\mathbf{X}^1(t=0) - \mathbf{X}^2(t=0)) \cdot \mathbf{n} \ge 0 \quad [6]
$$

les conditions d'impénétrabilité de Signorini sont alors données par :

$$
u_n + h \ge 0 , \quad r_n \ge 0 , \quad (u_n + h) r_n = 0 \tag{7}
$$
Ces conditions doivent être satisfaites à chaque instant $t \in \mathbf{I}$. Supposons maintenant que les solides soient initialement en contact sur une certaine partie de $\Gamma_c^\alpha$. Sur cette partie de $\Gamma_c^\alpha$, les conditions de Signorini deviennent alors :

$$ u_n \ge 0, \quad r_n \ge 0, \quad u_n r_n = 0 \qquad [8] $$

En général, à chaque instant $t \in \mathbf{I}$, la surface de contact potentiel $\Gamma_c^\alpha$ peut être découpée en deux parties disjointes : $+\Gamma_c^\alpha$ où les solides sont déjà en contact et $-\Gamma_c^\alpha$ où les solides ne sont pas en contact :

$$ \Gamma_c^\alpha = +\Gamma_c^\alpha \cup -\Gamma_c^\alpha \qquad [9] $$

Par rapport à $\Gamma_c^\alpha$, $+\Gamma_c^\alpha$ et $-\Gamma_c^\alpha$ changent au cours du temps $t$ et peuvent être vides à des instants $t \in \mathbf{I}$. Dans le cas de l'analyse dynamique, tel que les problèmes d'impact, les conditions de Signorini peuvent se formuler sur $+\Gamma_c^\alpha$, en termes de vitesse relative par :

$$ \dot{u}_n \ge 0, \quad r_n \ge 0, \quad \dot{u}_n r_n = 0 \quad \text{sur } +\Gamma_c^\alpha \qquad [10] $$

Lorsque $\dot{u}_n > 0$, les solides sont séparés alors qu'ils restent en contact lorsque $\dot{u}_n = 0$. La formulation [10] des conditions de Signorini peut être combinée avec les règles de glissement pour obtenir la loi complète de contact avec frottement applicable sur les parties en contact de $\Gamma_c^\alpha$. Cette loi complète exprime qu'il y a des vitesses possibles des solides qui satisfont la règle d'impénétrabilité, de non-adhesion et de glissement. Évidemment, pour un écart strictement positif ($u_n > 0$), la vitesse relative normale est arbitraire ($\dot{u}_n \in \mathbb{R}$) et la force normale de réaction est égale à zéro ($r_n = 0$). Les mouvements des corps qui ne sont pas en contact sont arbitraires jusqu'à ce que le contact soit établi. Ce choix est motivé par le fait que l'accent est mis sur la définition des évolutions admissibles des solides en contact où l'intégration temporelle doit être effectuée. Dans le reste du papier, un signe moins précède toujours la vitesse tangentielle relative $\dot{u}_t$ pour souligner sa direction opposée à la force de frottement.

Dans ce travail, on considère la loi classique de frottement sec de Coulomb. L'ensemble des forces admissibles, noté $K_\mu$, est défini par :

$$ K_\mu = \{ \mathbf{r} \in \mathbb{R}^3 \text{ tel que } |\mathbf{r}_t| - \mu r_n \le 0 \} \qquad [11] $$

où $\mu$ est le coefficient de frottement et $K_\mu$ le cône de Coulomb.

## 2.3. Loi complète de contact avec frottement

Sur la surface de contact $\Gamma_c^\alpha$, la règle de glissement peut être combinée avec les conditions de Signorini pour obtenir la loi complète de contact avec frottement qui indique les statuts possibles sur la zone de contact (adhérence, glissement, non-contact). Cette loi fortement non linéaire rend les problèmes de contact avec frotte-
ment parmi les plus difficiles à résoudre en mécanique. Une suite d'instructions du type "si...alors...sinon" peut être employée pour l'écrire analytiquement :

$$
\begin{array}{lll}
\text{si} & r_n = 0 \quad \text{alors} \quad \dot{u}_n > 0 & \text{non-contact} \\
\text{si} & \mathbf{r} \in K_{\mu}^{i} \quad \text{alors} \quad \dot{u}_n = 0 \quad \text{et} \quad \dot{\mathbf{u}}_t = \mathbf{0} & \text{adhérence} \\
\text{sinon} & (\mathbf{r} \in K_{\mu}^{f} \quad \text{et} \quad r_n > 0) & \\
& \left\{ \dot{u}_n = 0 \quad \text{et} \quad \exists \dot{\lambda} > 0 \quad \text{tel que} \quad -\dot{\mathbf{u}}_t = \dot{\lambda} \frac{\mathbf{r}_t}{|\mathbf{r}_t|} \right\} & \text{glissement} \hfill [12]
\end{array}
$$

où $K_{\mu}^{i}$ et $K_{\mu}^{f}$ représentent respectivement l'intérieur et la frontière de $K_{\mu}$. Le caractère
à valeurs multiples de la loi se situe dans la première et deuxième partie de la relation.
Si $r_n$ est nul alors $\dot{\mathbf{u}}$ est arbitraire, mais sa composante normale $\dot{u}_n$ devrait être positive.
En d'autres mots, un seul élément de $\mathbb{R}^3$ ($\mathbf{r} = \mathbf{0}$) est associé à un nombre infini
de vecteurs vitesse $\dot{\mathbf{u}} \in \mathbb{R}^3$. Les mêmes arguments peuvent être développés pour la
deuxième partie de la relation. La loi inverse, c'est-à-dire la relation $\mathbf{r}(-\dot{\mathbf{u}})$, peut être
écrite comme suit :

$$
\begin{array}{lll}
\text{si} & \dot{u}_n > 0 \text{ alors } r_n = 0 & \text{non-contact} \\
\text{si} & \mathbf{\dot{u}} = \mathbf{0} \text{ alors } \mathbf{r} \in K_\mu & \text{adhérence} \\
\text{sinon} & (\mathbf{\dot{u}} \in \mathbf{T} - \{\mathbf{0}\}) & \\
& \left\{ \dot{u}_n = 0 \text{ et } \mathbf{r}_t = \mu r_n \frac{-\mathbf{\dot{u}}_t}{| -\mathbf{\dot{u}}_t |} \right\} & \text{glissement} \hfill [13]
\end{array}
$$

La forme complète de la loi de contact avec frottement implique trois statuts possibles,
qui sont : non-contact, contact avec adhérence et contact avec glissement. Seul le
dernier cas produit de la dissipation d'énergie.

**3. Méthode du bipotentiel**

De Saxcé et Feng [DES 98] ont montré que la loi de contact [12] est équivalente à
l'inclusion différentielle suivante :

$$
-(\dot{\mathbf{u}}_t + (\dot{u}_n + \mu| - \dot{\mathbf{u}}_t|) \mathbf{n}) \in \partial \bigcup_{K_\mu} \mathbf{r} \quad [14]
$$

où $\bigcup_{K_{\mu}} \mathbf{r}$ est une fonction indicatrice de l'ensemble convexe fermé $K_{\mu}$ :

$$
\bigcup_{K_{\mu}} (\mathbf{r}) = \begin{cases} 0 & \text{si } \mathbf{r} \in K_{\mu} \\ +\infty & \text{sinon} \end{cases} \qquad [15]
$$

Le bipotentiel de contact est défini par :

$$
b_c(-\dot{\mathbf{u}}, \mathbf{r}) = \bigcup_{\mathbb{R}_-} (-\dot{u}_n) + \bigcup_{K_\mu} (\mathbf{r}) + \mu r_n | -\dot{\mathbf{u}}_t| \quad [16]
$$
où $\mathbb{R}_- = ] - \infty, 0]$ est l'ensemble des nombres réels négatifs ou nuls.

Afin d'éviter les potentiels non différentiables qui apparaissent en mécanique non linéaire, tel que les problèmes de contact, il est d'usage d'employer la méthode du Lagrangien augmenté [ALA 91, SIM 92, DES 91, KLA 92]. Pour le bipotentiel de contact $b_c$, donné par [16], à condition que $\dot{u}_n \ge 0$ et $\mathbf{r} \in K_\mu$, nous avons :

$$ \forall \mathbf{r}' \in K_\mu, \quad \varrho\mu(r'_n - r_n) - |\dot{\mathbf{u}}_t| + (\mathbf{r} - (\mathbf{r} - \varrho\dot{\mathbf{u}}))(\mathbf{r}' - \mathbf{r}) \ge 0 \quad [17] $$

où $\varrho$ est un paramètre intervenant dans la solution. Dans nos calculs, $\varrho$ est pris comme étant égal au plus grand terme diagonal de matrice de rigidité de contact local. Tenant compte de la décomposition [3], l'inégalité suivante doit être satisfaite :

$$ \mathbf{r}' \in K_\mu, \quad (\mathbf{r} - \boldsymbol{\tau}) \cdot (\mathbf{r}' - \mathbf{r}) \ge 0 \quad [18] $$

où les réactions de contact augmentées $\boldsymbol{\tau}$ sont définies par :

$$ \boldsymbol{\tau} = \mathbf{r} - \varrho(\dot{\mathbf{u}}_t + (\dot{u}_n + \mu| - \dot{\mathbf{u}}_t|)\mathbf{n}) \quad [19] $$

L'inégalité [18] signifie que $\mathbf{r}$ est la projection de $\boldsymbol{\tau}$ sur le cône de Coulomb :

$$ \mathbf{r} = \mathrm{proj}(\boldsymbol{\tau}, K_\mu) \quad [20] $$

pour la solution numérique de l'équation implicite [20], l'algorithme d'Uzawa est employé, celui-ci conduit à un processus itératif impliquant une étape de prédiction/correction :

$$ \begin{array}{l} \text{Prédiction} \quad \boldsymbol{\tau}^{i+1} = \mathbf{r}^i - \varrho^i (\dot{\mathbf{u}}_t^i + (\dot{u}_n^i + \mu| - \dot{\mathbf{u}}_t^i|) \mathbf{n}) \\ \text{Correction} \quad \mathbf{r}^{i+1} = \mathrm{proj}(\boldsymbol{\tau}^{i+1}, K_\mu) \end{array} \quad [21] $$

Il est important de noter que, dans cet algorithme, le contact unilatéral et le frottement sont couplés par l'intermédiaire du bipotentiel. Un autre point essentiel de la méthode du bipotentiel est que le correcteur peut être analytiquement trouvé en ce qui concerne les trois statuts possibles du contact (figure 2) : $\boldsymbol{\tau} \subset K_\mu$ (contact avec adhérence), $\boldsymbol{\tau} \subset K_\mu^*$ (pas de contact) and $\boldsymbol{\tau} \subset \mathbb{R}^3 - K_\mu \cup K_\mu^*$ (contact avec glissement). $K_\mu^*$ est le cône dual de $K_\mu$. Cette étape de correction est explicitement donnée comme suit :

$$ \text{si} \quad \mu|\boldsymbol{\tau}_t^{i+1}| < -\tau_n^{i+1} \quad \text{alors} \quad \mathbf{r}^{i+1} = 0 \quad \text{non-contact} $$

$$ \text{si} \quad |\boldsymbol{\tau}_t^{i+1}| < \mu \tau_n^{i+1} \quad \text{alors} \quad \mathbf{r}^{i+1} = \boldsymbol{\tau}^{i+1} \quad \text{adhérence} \quad [22] $$

$$ \text{sinon} \quad \mathbf{r}^{i+1} = \boldsymbol{\tau}^{i+1} - \frac{(|\boldsymbol{\tau}_t^{i+1}| - \mu \tau_n^{i+1})}{(1 + \mu^2)} \left( \frac{\boldsymbol{\tau}_t^{i+1}}{|\boldsymbol{\tau}_t^{i+1}|} + \mu \mathbf{n} \right) \quad \text{glissement} $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Correction sur le cône de Coulomb selon les statuts du contact</figcaption></figure>

Il est important de souligner que cette formule explicite est valable à la fois pour les problèmes de contact 2D et 3D avec frottement de Coulomb et permet d'obtenir des résultats très stables et précis.

# 4. Formulation éléments finis

## 4.1. Formulation du Lagrangien total

Dans l'analyse linéaire, on suppose une relation linéaire entre les déformations et les déplacements. Cependant, s'il y a de grands déplacements et de grandes rotations, comme dans le cas des problèmes de contact en dynamique, la relation non linéaire entre les déformations et les déplacements ne peut pas être ignorée. L'analyse non linéaire géométrique peut être décrite en employant la formulation du Lagrangien total ou du Lagrangien actualisé. La formulation du Lagrangien total construit la matrice de rigidité tangente par rapport à la configuration initiale, alors que la formulation du Lagrangien actualisé construit la matrice de rigidité tangente à partir de la configuration courante. Dans la formulation du Lagrangien total, la configuration initiale demeure constante. Ceci simplifie le calcul [CRI 91, SIM 98]. Par conséquent, la formulation du Lagrangien total a été choisie dans ce travail pour la discrétisation par éléments finis. Afin de décrire la non-linéarité géométrique, le tenseur du gradient des déformations est défini par :

$$ \Phi = \mathrm{Id} + \nabla u \quad [23] $$

où $\mathrm{Id}$ est le tenseur unité et $\nabla u$ le tenseur gradient des déplacements. En raison des grands déplacements et des grandes rotations, le tenseur de Green-Lagrange est utilisé pour exprimer la relation non linéaire entre les déformations et les déplacements. Nous
notons C le tenseur des dilatations ou tenseur Cauchy-Green droit ($C = \Phi^T \Phi$). Le tenseur des déformations de Green-Lagrange E est alors défini par :

$$E = \frac{1}{2}(C - I) \quad [24]$$

Dans le cadre de la méthode des éléments finis et d'après les équations [23] et [24], le tenseur de Green-Lagrange inclut formellement les termes linéaires et non linéaires en fonction des déplacement nodaux :

$$E = (B_L + \frac{1}{2}B_{NL}(u))u \quad [25]$$

où $B_L$ et $B_{NL}(u)$ sont les matrices qui relient respectivement les termes linéaires et non linéaires des déformations aux déplacements nodaux. D'après [25], la forme incrémentale de la relation déformation-déplacement est :

$$\delta E = (B_L + B_{NL}(u))\delta u \quad [26]$$

Dans le cas des lois élastiques ou hyper-élastiques, il existe un potentiel élastique $W$ (ou densité d'énergie) qui est une fonction scalaire du tenseur des déformations, dont la dérivée par rapport aux déformations correspond aux contraintes. Ce qui peut s'exprimer par :

$$S = \frac{\partial W}{\partial E} = 2 \frac{\partial W}{\partial C} \quad [27]$$

où $S$ est le second tenseur des contraintes de Piola-Kirchhoff. Dans le cas particulier des modèles isotropes de Saint-Venant-Kirchhoff, nous avons :

$$W = \frac{1}{2} E : DE \quad [28]$$

Alors $S$ peut s'écrire :

$$S = DE \quad [29]$$

où $D$ est le tenseur d'élasticité.

Dans un cadre général, le travail virtuel $\delta U$ est donné par :

$$\delta U = M \dot{u} \delta u + A \dot{u} \delta u + \int_{V_0} S \delta E dV - F_{ext} \delta u - R \delta u = 0 \quad [30]$$

où $V_0$ est le volume de la configuration initiale, $F_{ext}$ le vecteur des chargements extérieurs, $R$ le vecteur des réactions de contact, $M$ la matrice de masse, $A$ la matrice d'amortissement, $\dot{u}$ le vecteur vitesse et $\ddot{u}$ le vecteur accélération. En remplaçant $\delta E$ donné par [26] dans [30], on obtient :

$$\delta U = M \ddot{u} \delta u + A \dot{u} \delta u + \int_{V_0} (B_L + B_{NL})^T S \delta u dV - F_{ext} \delta u - R \delta u = 0 \quad [31]$$
Le vecteur des forces internes est alors défini par :

$$ \mathbf{F}_{int} = \int_{V_0} (\mathbf{B}_L + \mathbf{B}_{NL}(\mathbf{u}))^T \mathbf{S} \, dV \quad [32] $$

Puisque $\delta \mathbf{u}$ est arbitraire, on obtient le système d'équations non linéaires suivant :

$$ \mathbf{M}\ddot{\mathbf{u}} + \mathbf{A}\dot{\mathbf{u}} + \mathbf{F}_{int} - \mathbf{F}_{ext} - \mathbf{R} = 0 \quad [33] $$

Il est à noter que la rigidité est prise en considération par le vecteur des forces internes $\mathbf{F}_{int}$. L'équation [33] peut être réécrite sous la forme :

$$ \mathbf{M}\ddot{\mathbf{u}} = \mathbf{F} + \mathbf{R} \quad \text{avec} \quad \mathbf{F} = \mathbf{F}_{ext} - \mathbf{F}_{int} - \mathbf{A}\dot{\mathbf{u}} \quad [34] $$

avec les conditions initiales à $t = 0$ :

$$ \dot{\mathbf{u}} = \dot{\mathbf{u}}_0 \quad \text{et} \quad \mathbf{u} = \mathbf{u}_0 \quad [35] $$

La dérivée de $\mathbf{F}_{int}$ par rapport aux déplacements nodaux $\mathbf{u}$ donne la matrice de rigidité tangente :

$$ \mathbf{K} = \frac{\partial \mathbf{F}_{int}}{\partial \mathbf{u}} = \int_{V_0} \left( \frac{\partial \mathbf{S}}{\partial \mathbf{u}} (\mathbf{B}_L + \mathbf{B}_{NL}(\mathbf{u})) + \mathbf{S} \frac{\partial \mathbf{B}_{NL}(\mathbf{u})}{\partial \mathbf{u}} \right) dV \quad [36] $$

De plus, en utilisant [26] et [29], la matrice de rigidité tangente peut se décomposer en trois matrices qui sont la matrice de rigidité élastique $\mathbf{K}_e$, la matrice de rigidité géométrique (ou la matrice de rigidité des contraintes initiales) $\mathbf{K}_\sigma$ et la matrice de rigidité des déplacements initiaux $\mathbf{K}_u$ :

$$ \mathbf{K} = \mathbf{K}_e + \mathbf{K}_\sigma + \mathbf{K}_u \quad [37] $$

où

$$ \mathbf{K}_e = \int_{V_0} \mathbf{B}_L^T \mathbf{D} \mathbf{B}_L \, dV \quad [38] $$

$$ \mathbf{K}_\sigma = \int_{V_0} \mathbf{S} \frac{\partial \mathbf{B}_{NL}}{\partial \mathbf{u}} \, dV \quad [39] $$

$$ \mathbf{K}_u = \int_{V_0} (\mathbf{B}_L^T \mathbf{D} \mathbf{B}_{NL} + \mathbf{B}_{NL}^T \mathbf{D} \mathbf{B}_L + \mathbf{B}_{NL}^T \mathbf{D} \mathbf{B}_{NL}) \, dV \quad [40] $$

## 4.2. Algorithme d'intégration du premier ordre

Nous pouvons maintenant intégrer l'équation [34] entre deux instants consécutifs $t$ et $t + \Delta t$. Pour réaliser cette opération, il est classique d'utiliser l'algorithme
du second ordre de Newmark. Cependant, dans les problèmes d'impact, les approximations d'ordre supérieur ne signifient pas nécessairement une meilleure précision. En effet, lors du changement brutal des statuts de contact, la vitesse et l'accélération ne sont pas continues, et la régularité excessive peut conduire à des erreurs importantes [ARM 98]. Pour cette raison, Jean [JEA 89, JEA 99] a proposé un algorithme du premier ordre qui est employé dans ce travail. Cet algorithme est basé sur les approximations suivantes :

$$ \mathbf{u}^{t+\Delta t} - \mathbf{u}^t = \Delta t \left( (1-\theta) \dot{\mathbf{u}}^t + \theta \dot{\mathbf{u}}^{t+\Delta t} \right) \quad [41] $$

$$ \int_t^{t+\Delta t} \mathbf{M} \, d\dot{\mathbf{u}} = \mathbf{M} (\dot{\mathbf{u}}^{t+\Delta t} - \dot{\mathbf{u}}^t) \quad [42] $$

$$ \int_t^{t+\Delta t} \mathbf{F} \, dt = \Delta t \left( (1-\xi) \mathbf{F}^t + \xi \mathbf{F}^{t+\Delta t} \right) \quad [43] $$

$$ \int_t^{t+\Delta t} \mathbf{R} \, dt = \Delta t \mathbf{R}^{t+\Delta t} \quad [44] $$

où $0 \le \xi \le 1$; $0 \le \theta \le 1$. Pour le processus itératif, les valeurs au temps $t + \Delta t$ sont remplacées par les valeurs à l'itération $i + 1$; soit par exemple, $\mathbf{F}^{t+\Delta t} = \mathbf{F}^{i+1}$. L'approximation standard de $\mathbf{F}^{i+1}$ donne :

$$ \mathbf{F}^{i+1} = \mathbf{F}_{\text{int}}^i + \frac{\partial \mathbf{F}}{\partial \mathbf{u}} (\mathbf{u}^{i+1} - \mathbf{u}^i) + \frac{\partial \mathbf{F}}{\partial \dot{\mathbf{u}}} (\dot{\mathbf{u}}^{i+1} - \dot{\mathbf{u}}^i) = \mathbf{F}_{\text{int}}^i - \mathbf{K}^i \Delta \mathbf{u} - \mathbf{A}^i \Delta \dot{\mathbf{u}} \quad [45] $$

Ainsi, on obtient la forme recursive en termes de déplacement :

$$ \begin{aligned} \bar{\mathbf{K}}^i \Delta \mathbf{u} &= \bar{\mathbf{F}}^i + \bar{\mathbf{F}}_{\text{acc}}^i + \mathbf{R}^{i+1} \\ \mathbf{u}^{i+1} &= \mathbf{u}^i + \Delta \mathbf{u} \end{aligned} \quad [46] $$

où les différents quantités effectives sont données par :

$$ \bar{\mathbf{K}}^i = \xi \mathbf{K}^i + \frac{\xi}{\theta \Delta t} \mathbf{A}^i + \frac{1}{\theta \Delta t^2} \mathbf{M}^i \quad [47] $$

$$ \bar{\mathbf{F}}_{\text{acc}}^i = -\frac{1}{\theta \Delta t^2} \mathbf{M}^i \{\mathbf{u}^i - \mathbf{u}^t - \Delta t \dot{\mathbf{u}}^t\} \quad [48] $$

$$ \bar{\mathbf{F}}^i = (1 - \xi) (\mathbf{F}_{\text{int}}^t + \mathbf{F}_{\text{ext}}^t) + \xi (\mathbf{F}_{\text{int}}^i + \mathbf{F}_{\text{ext}}^{t+\Delta t}) \quad [49] $$

À la fin de chaque pas de temps, la vitesse est réactualisée par :

$$ \dot{\mathbf{u}}^{t+\Delta t} = \left(1 - \frac{1}{\theta}\right) \dot{\mathbf{u}}^t + \frac{1}{\theta \Delta t} (\mathbf{u}^{t+\Delta t} - \mathbf{u}^t) \quad [50] $$
En prenant $\theta = \frac{1}{2}$, ce schéma correspond à la règle du trapèze implicite, qui est également analogue à la méthode développée par Tamma et Namburu [TAM 90] dans laquelle l'accélération n'a pas besoin d'être calculée. Simo et Wong [SIM 91] ont montré que ce schéma préserve l'énergie totale et l'équilibre pour les problèmes dynamiques sans contact. Les exemples vont montrer que c'est également le cas avec les problèmes d'impact sans frottement.

Il faut noter que l'équation [46] est fortement non linéaire, en raison des grandes rotations et des grands déplacements du solide, comme c'est le cas pour les problèmes de contact et d'impact de système multicorps. En outre, comme il est précisé plus haut, la loi de contact avec frottement est habituellement décrite par une inégalité et le potentiel de contact est même non différentiable. Afin de résoudre cette équation en tenant compte de toutes les non-linéarités en même temps, Feng [FEN 95] a proposé une stratégie qui consiste à traiter séparément les non-linéarités afin de réduire la complexité du calcul et d'améliorer la stabilité numérique. Comme $\Delta u$ et $\mathbf{R}$ sont inconnus, l'équation [46] ne peut être résolue directement. Dans un premier temps, le vecteur $\mathbf{R}$ est calculé à l'aide de la méthode du bipotentiel dans un système réduit qui ne concerne que les nœuds de contact. Ensuite, le vecteur $\Delta u$ peut être déterminé en considérant les réactions de contact comme un chargement extérieur. Il est important de souligner que contrairement à la méthode de pénalisation ou la méthode des multiplicateurs de Lagrange, la méthode du bipotentiel ne modifie pas la matrice de rigidité globale, ni n'augmente le nombre de degrés de liberté. La conséquence de cette intéressante propriété est la facilité à implanter cette méthode de résolution pour les problèmes de contact avec frottement au sein d'un code de calcul par éléments finis. De plus, en raison de la séparation des non-linéarités, la méthode est plus stable et plus efficace.

### 4.3. Calcul d'énergie

Après détermination du champ de déplacement et de vitesse, nous pouvons calculer différentes énergies. Toute l'énergie élastique des solides (discrétisés par $n_{el}$ éléments finis) est obtenue par la relation suivante :

$$E_e = \sum_{e=1}^{n_{el}} \int_{\Omega_e} W_e \, d\Omega \quad [51]$$

L'énergie cinétique peut être calculée au niveau global par :

$$E_k = \frac{1}{2} \dot{\mathbf{u}}^T \mathbf{M} \dot{\mathbf{u}} \quad [52]$$

Finalement, l'énergie totale du système de solides correspond à :

$$E_t = E_e + E_k \quad [53]$$

Les cas traités dans la suite sont des problèmes homogènes de Neumann, caractérisés par aucun déplacement imposé sur la frontière et aucun chargement extérieur. En
outre, si le contact sans frottement est considéré, l'énergie totale devrait être conser-
vée. L'algorithme proposé préserve parfaitement cette propriété fondamentale de con-
servation de l'énergie comme le montrent les exemples numériques.

**5. Résultats numériques**

Les algorithmes développés ont été implantés et testés dans le code de calcul par
éléments finis FER/Impact [FEN 00]. Par ailleurs, des exemples d'application, dans
des cas statiques ou quasi statiques, ont été effectués [FEN 95, FEN 98, FEN 03].

Deux exemples vont être décrits, dans cet article, afin d'illustrer les performances
de l'algorithme pour la résolution de problèmes de contact/impact avec frottement en
dynamique. Pour tous les cas, on suppose qu'il n'y a pas d'amortissement à l'exception du frottement de Coulomb entre les surfaces de contact, c'est-à-dire **A** = 0 dans
les équations [34] et [47].

Le premier exemple, déjà étudié par Kwak *et al.* [KO 92, KIM 96] en utilisant une méthode de programmation mathématique, va servir à valider l’algorithme tout en démontrant son efficacité. Une plaque vient impacter suivant une direction oblique une surface rigide puis elle rebondit. La plaque est discrétisée par des éléments quadrilatéraux en contrainte plane (figure 3). Deux maillages sont utilisés pour l’étude : M1 (54 nœuds et 37 éléments) et M2 (292 nœuds et 258 éléments). Afin d’effectuer des comparaisons, le maillage M1 est identique à celui utilisé dans [KIM 96].

<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Impact oblique d'une plaque élastique : géométrie et maillages*</figcaption></figure>

Les caractéristiques de cet exemple sont : Module d'Young E = 10⁷ Pa, coefficient de Poisson v = 0,25, masse volumique ρ = 1000 kg/m³, coefficient de frottement μ = 0,1, vitesse initiale : vₓ = 3 m/s, vᵧ = -5 m/s. Les dimensions de la plaque sont : L = 0,04 m, H = 0,08 m, rayon R = 0,101 m, épaisseur e = 0,01 m.
Le temps de simulation est $3 \times 10^{-3}$ s et les paramètres du schéma sont : $\Delta t = 10^{-5}$ s, $\xi = \theta = 0,5$.

Dans un premier temps, on considère le maillage M1. La figure 4 montre clairement que l’énergie totale est dissipée par frottement, soulignons au passage que cette énergie dissipée est calculée de manière quantitative. Par ailleurs, il est peut-être intéressant de savoir si l’énergie dissipée est proportionnelle au coefficient de frottement. La simulation numérique montre que ce n’est pas le cas comme on peut le constater sur la figure 5. Ce résultat peut s’interpréter simplement. En effet, lorsque le coefficient de frottement augmente, les forces de frottement augmentent également, mais le glissement tangentiel va quant à lui diminuer. La dissipation d’énergie dépendant à la fois des forces de frottement et du glissement tangentiel au niveau des nœuds de contact traduit le comportement observé sur la figure 5.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Évolution de l’énergie au cours du temps</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Variation de l’énergie dissipée en fonction du coefficient de frottement</figcaption></figure>
À titre indicatif, les temps CPU de la méthode proposée et ceux obtenus par Kwak *et al.* sont reportés dans le tableau 1.

<table><caption>**Tableau 1.** Comparaison des temps CPU</caption><thead><tr><td>Méthode</td><td>Machine</td><td>Temps CPU (s)</td></tr></thead><tbody><tr><td>Ko & Kwak [KO 92]</td><td>CRAY 2S/4-128</td><td>19 000</td></tr><tr><td>Kim & Kwak [KIM 96]</td><td>HP 720</td><td>430</td></tr><tr><td>Méthode proposée</td><td>PC Pentium 4/2.8 GHz</td><td>7</td></tr></tbody></table>

On analyse maintenant l'influence du pas de temps sur les résultats. Pour cela, on considère les pas de temps suivants : $\Delta t = 10^{-3}$ s (A), $\Delta t = 10^{-4}$ s (B), $\Delta t = 10^{-5}$ s (C) et $\Delta t = 10^{-6}$ s (D). Les composantes ($U_x, U_y$) de la trajectoire du point P (figure 3) sont tracées en fonction des différents pas de temps sur la figure 6. Globalement, on ne constate pas de différences importantes. Évidemment, lorsque le pas de temps devient grand, on commence à percevoir des différences sur la trajectoire.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Influence du pas de temps sur les déplacements du point P</figcaption></figure>

Par ailleurs, pour assurer la stabilité inconditionnelle du schéma d'intégration, le paramètre $\xi$ doit être supérieur ou égal à 0,5. Lorsque l'on augmente la valeur de ce paramètre, on constate une diminution d'énergie totale par dissipation numérique (figure 7).
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Influence du coefficient ξ sur l'énergie totale</figcaption></figure>

Enfin, on s'intéresse à l'influence de la finesse du maillage. Les figures 8 et 9 représentent, pour le point P, respectivement les déplacements et les contraintes en fonction du maillage (M1 et M2).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Influence du maillage sur les déplacements du point P</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Influence du maillage sur les contraintes au point P</figcaption></figure>

On propose un deuxième exemple test (sans unité) qui simule l'impact d'un disque entre deux plans rigides espacés de 6 (figure 10). Le disque a un diamètre de 2 et une épaisseur de 1. Les propriétés du matériau (modèle de Saint-Venant-Kirchhoff) sont : Module d'Young $E = 1000$, coefficient de Poisson $v = 0,45$, et masse volumique $\rho = 1$. La vitesse initiale du disque est : $v_x = 2, v_y = -2$. Le temps de simulation est 8 et les paramètres du schéma sont : $\Delta t = 10^{-3}$, $\xi = \theta = 0,5$. Cet exemple montre clairement l'influence du frottement sur la cinématique du disque et sur l'évolution de l'énergie. Les figures (10, 11) représentent l'évolution du disque à différents instants ($t = 0, t = 0,1588, t = 3,888, t = 6,258$ et $t = 7,498$) respectivement dans le cas du contact sans frottement et avec frottement ($\mu = 0,5$). Dans le premier cas, l'évolution est prévisible, le disque se déforme sous l'impact et se déplace sans rotation. Par contre dans le second cas, l'évolution du disque est moins prévisible, puisque le disque revient vers sa position initiale en raison de sa rotation induite par le frottement.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Rebonds sans frottement d'un disque entre deux plans rigides</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Rebonds avec frottement d'un disque entre deux plans rigides</figcaption></figure>

La figure 12 représente la distribution des contraintes de Von-Mises dans le disque au moment du premier impact dans les deux cas. Dans le cas sans frottement (figure de gauche), la distribution est comme prévue symétrique, mais elle ne l'est pas du tout dans le cas avec frottement (figure de droite).

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Distribution des contraintes de Von-Mises lors d'un impact sans frottement et avec frottement</figcaption></figure>

Les figures 13 et 14 représentent respectivement les énergies dans le cas du contact sans frottement et avec frottement. Dans le cas sans frottement, on constate que l'énergie totale $E_t$ est parfaitement conservée. Pour le cas avec frottement, l'énergie totale diminue après chaque impact en raison des frottements.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** *Évolution de l'énergie, cas sans frottement*</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** *Évolution de l'énergie, cas avec frottement*</figcaption></figure>

**6. Conclusion**

Dans cet article, nous avons présenté l'extension récente de la méthode du bipotentiel à l'analyse dynamique des problèmes bidimensionnels de contact avec frottement de Coulomb. L'algorithme développé a été décrit et étudié numériquement pour deux problèmes en utilisant différents coefficients de frottement. Les résultats numériques montrent que :
– L’algorithme développé préserve quasi parfaitement la propriété de conservation d’énergie du contact sans frottement.

– L’algorithme développé permet de déterminer quantitativement la dissipation physique d’énergie dus au frottement.

– L’énergie dissipée n’est pas proportionnelle au coefficient de frottement.

– L’algorithme développé est simple et efficace :

- aucune modification de la matrice globale de rigidité ;

- aucune régularisation des lois de contact et de frottement ;

- calcul précis des forces de contact dans un système réduit ;

- intégration du premier ordre au lieu de l’intégration de second ordre ou d’ordre supérieur.

Il est convient de noter que cet algorithme permet également de traiter les problèmes de contact/impact entre plusieurs corps déformables [FEN 03, FEN 04]. De plus, nous pensons que cette approche pourrait facilement être étendue aux problèmes dynamiques tridimensionnels de contact comprenant des lois constitutives non linéaires et des modèles de frottement plus complexes [HJI 04]. Ce travail est en cours de développement.

## 7. Bibliographie

[ALA 91] ALART P., CURNIER A., « A mixed formulation for frictional contact problems prone to Newton like solution methods », *Comp. Meth. Appl. Mech. Engng.*, vol. 92, 1991, p. 353-375.

[ARM 98] ARMERO F., PETOCZ E., « Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems », *Comp. Meth. Appl. Mech. Engng.*, vol. 158, 1998, p. 269-300.

[CHA 86] CHAUDHARY A. B., BATHE K. J., « A solution method for static and dynamic analysis of three dimensional contact problems with friction », *Computers and Structures*, vol. 24, 1986, p. 855-873.

[CHR 98] CHRITENSEN P. W., KLARBRING A., PANG J. S., STRÖMBERG N. N., « Formulation and comparison of algorithms for frictional contact problems », *Int. J. Num. Meth. Engng.*, vol. 42, 1998, p. 145-173.

[CRI 91] CRISFIELD M. A., *Non-linear finite element analysis of solid and structures*, Wiley, 1991.

[DES 91] DE SAXCÉ G., FENG Z.-Q., « New inequality and functional for contact with friction : The implicit standard material approach », *Mech. Struct. & Mach.*, vol. 19, 1991, p. 301-325.

[DES 98] DE SAXCÉ G., FENG Z.-Q., « The bi-potential method : a constructive approach to design the complete contact law with friction and improved numerical algorithms », *Mathematical and Computer Modeling*, vol. 28(4-8),, 1998, p. 225-245, Special issue : Recent Advances in Contact Mechanics.
[FEN 95] FENG Z.-Q., « 2D or 3D frictional contact algorithms and applications in a large deformation context », *Comm. Numer. Meth. Engng.*, vol. 11, 1995, p. 409-416.

[FEN 98] FENG Z.-Q., « Some test examples of 2D and 3D contact problems involving Coulomb friction and large slip », *Mathematical and Computer Modeling*, vol. 28(4-8), 1998, p. 469-477, Special issue : Recent Advances in Contact Mechanics.

[FEN 00] FENG Z.-Q., « http://gmfe16.cemif.univ-evry.fr:8080/~feng/FerSystem.html », 2000.

[FEN 03] FENG Z.-Q., PEYRAUT F., LABED N., « Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies », *Int. J. Engn. Science*, vol. 41, 2003, p. 2213-2225.

[FEN 04] FENG Z.-Q., MAGNAIN B., CROS J.-M., JOLI P., « Energy dissipation by friction in dynamic multibody contact problems », YAO Z.-H., YUAN M.-W., ZHONG W.-X., Eds., *Computational Mechanics*, Beijing, China, sep 2004, WCCM VI in conjunction with APCOM04, Springer.

[HJI 04] HJIAJ M., FENG Z.-Q., DE SAXCÉ G., MRÓZ Z., « There-dimensional finite element computations for frictional contact problems with non-associated sliding rule », *Int. J. Num. Meth. Engng.*, vol. 60, 2004, p. 2045-2076.

[HUG 76] HUGHES T. J. R., TAYLOR R. L., SACKMAN J. L., CURNIER A., KANOKNUKUL-CHAI W., « A finite element method for a class of contact-impact problems », *Comp. Meth. Appl. Mech. Engng.*, vol. 8, 1976, p. 249-276.

[JEA 89] JEAN M., « Dynamics with partially elastic shocks and dry friction : double scale method and numerical approach », *4th Meeting on unilateral problems in structural analysis*, 1989, Capri.

[JEA 99] JEAN M., « The non-smooth contact dynamics method », *Comp. Meth. Appl. Mech. Engng.*, vol. 177, 1999, p. 235-257.

[KIK 88] KIKUCHI N., ODEN J. T., *Contact problems in elasticity : A study of variational inequalities and finite elements*, Philadelphia : SIAM, 1988.

[KIM 96] KIM J. O., KWAK B., « Dynamic analysis of two-dimensional frictional contact by linear complementarity problem formulation », *Int. J. Solids and Structures*, vol. 33, 1996, p. 4605-4624.

[KLA 92] KLARBRING A., « Mathematical programming and augmented Lagrangian methods for frictional contact problems », CURNIER A., Ed., *Contact Mechanics Int. Symp.*, 1992, PPUR.

[KLA 93] KLARBRING A., « Mathematical programming in contact problems », ALIABDALI M., BREBBIA C., Eds., *Computational methods in contact mechanics*, Southampton : Computational Mechanics Publications, 1993, p. 233-263.

[KO 92] KO S. H., KWAK B., « Frictional dynamic contact analysis using finite element nodal displacement description », *Computers & Structures*, vol. 42, 1992, p. 797-807.

[LAU 97] LAURSEN T., CHAWLA V., « Design of energy conserving algorithms for frictionless dynamic contact problems », *Int. J. Num. Meth. Engn.*, vol. 40, 1997, p. 863-886.

[LAU 02] LAURSEN T. A., LOVE G. R., « Improved implicit integrators for transient impact problems geometric admissibility within the conserving framework », *Int. J. Num. Meth. Engng.*, vol. 53, 2002, p. 245-274.
[PAR 89] PARISCH H. A., « Consistent tangent stiffness matrix for three-dimensional non-linear contact analysis », *Int. J. Num. Meth. Engng.*, vol. 28, 1989, p. 1803-1812.

[SIM 91] SIMO J. C., WONG K. K., « Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum », *Int. J. Numer. Meth. Engng.*, vol. 31, 1991, p. 19-52.

[SIM 92] SIMO J. C., LAURSEN T. A., « An augmented Lagrangian treatment of contact problems involving friction », *Computers & Structures*, vol. 42, 1992, p. 97-116.

[SIM 98] SIMO J. C., HUGHES T. J. R., *Computational inelasticity*, Springer-Verlag, New York, 1998.

[TAM 90] TAMMA K. K., NAMBURU R. R., « A robust self - starting explicit computational methodology for structural dynamic applications : architecture and representations », *Int. J. Numer. Meth. Engng.*, vol. 30, 1990, p. 1441-1454.

[WRI 90] WRIGGERS P., VU VAN T., STEIN E., « Finite element formulation of large deformation impact contact problems with friction », *Computers & Structures*, vol. 37, 1990, p. 319-331.

[WRI 95] WRIGGERS P., « Finite element algorithms for contact problems », *Arch. Comput. Meth. Engng.*, vol. 2, 1995, p. 1-49.

[WRI 02] WRIGGERS P., *Computational contact mechanics*, John Wiley & Sons, 2002.

[ZHO 93] ZHONG Z. H., *Finite element procedures in contact-impact problems*, Oxford University Press, 1993.

[ZIE 84] ZIENKIEWICZ O. C., WOOD W. L., HINE L. W., TAYLOR R. L., « A unified set of single step algorithms. Part 1. General formulation and application », *Int. J. Num. Meth. Engng.*, vol. 20, 1984, p. 1529-1552.