# Méthodes duales pour le contact frottant

François Kuss, Frédéric Lebon

► To cite this version:

François Kuss, Frédéric Lebon. Méthodes duales pour le contact frottant. XVIIe Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-00108990

HAL Id: hal-00108990

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

Submitted on 29 Jan 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
# Méthodes duales pour le contact frottant

François Kuss et Frédéric Lebon

Laboratoire de Mécanique et d'Acoustique
Université de Provence Aix-Marseille 1
31 Chemin Joseph Aiguier 13402 Marseille Cedex 20
(kuss,lebon)@lma.cnrs-mrs.fr

## Résumé :

On présente deux algorithmes destinés à la résolution de problèmes de contact frottant. Ces algorithmes sont basés sur la formulation duale ou complémentaire du problème de frottement. Des résultats numériques sont présentés et analysés.

## Mots-clefs :

## Mécanique du contact ; éléments finis ; énergie complémentaire

### 1 Introduction

La mécanique du contact intervient dans de très nombreux problèmes d'ingénierie (automobile, aéronautique, énergie, génie civil, biomécanique, géophysique, ...) et soulève de nombreux problèmes théoriques et numériques. La résolution de ces problèmes dans des logiciels industriels est encore mal maîtrisée surtout s'ils sont couplés avec d'autres non linéarités. Généralement, les méthodes de résolutions sont basées sur des formulations en déplacements (ou en vitesses) ou mixtes déplacements/forces de contact. La stratégie duale ou complémentaire faisant intervenir une formulation purement en forces de contact a par contre été fort peu développée jusqu'à présent [1, 2, 3] malgré des résultats prometteurs. Cette stratégie se heurte à deux écueils : le choix des éléments finis vérifiant les équations d'équilibre et les méthodes itératives de résolution des systèmes non linéaires obtenus. Dans cet article, on propose d'utiliser une famille d'éléments faisant intervenir les fonctions d'Airy et des algorithmes de résolutions basés sur des notions de relaxation ponctuelle ou par blocs. Des résultats numériques sont donnés et des comparaisons avec des méthodes traditionnelles sont présentées.

### 2 Quelques rappels sur les éléments finis en contraintes

#### 2.1 Formulation duale du problème de contact avec frottement

On s'intéresse à l'équilibre d'un corps élastique $\Omega$, soumis à des déplacements imposés, à des forces de contact, volumiques et surfaciques respectivement sur $\Gamma_D$, $\Omega$, $\Gamma_C$ et $\Gamma_F$. La formulation duale du problème d'élasticité en présence de contact et de frottement est donnée par :

$$ P_c \begin{cases} \text{Trouver } \sigma \in K_{C(-\sigma_N)} \text{ telle que} \\ b(\sigma, \tau - \sigma) \ge l(\tau - \sigma) \quad \forall \tau \in K_{C(-\sigma_N)} \end{cases} \quad (1) $$
Avec :

$$b(\sigma, \tau) = \int_{\Omega} A_{ijkl}\sigma_{ij}\tau_{kl}dx$$

$$l(\tau) = \int_{\Gamma_D} \tau_{ij} n_j u_{0i} ds$$

et

$$H = \{ \tau \in L^2(\Omega, M); \tau_{ij,j} + F_i = 0 \text{ dans } \Omega \text{ et } \tau n = f \text{ sur } \Gamma_F \} \quad (2)$$

$$K_{C(\eta)} = \{ \tau \in H; \tau_N \le 0 \text{ sur } \Gamma_C; \|\tau_T\| \le \mu\eta \text{ sur } \Gamma_C \}$$

Avec M l'espace des tenseurs symétriques du second ordre et $K_{C(-\sigma_N)}$, l'ensemble convexe des contraintes admissibles. Cet ensemble dépendant de la contrainte elle même, le problème formulé en contraintes ne peut alors pas être réduit à un problème de minimisation.

## 2.2 Discrétisation du champ de contraintes

Le champ de contraintes sur un élément doit satisfaire les conditions d'équilibre de volume et de surfaces : $\tau_{ij,j} + F_i = 0$ dans $\Omega$ et $\tau n = f$ sur $\Gamma_F$. De manière à satisfaire la condition d'équilibre de volume, on définit en deux dimensions le champ de contraintes comme dérivant d'une fonction de Airy ($A(x,y)$), telle que, dans un repère cartésien orthonormé :

$$\tau_{\alpha\beta} = e_{\alpha\gamma}e_{\beta\delta}A(x,y)_{,\gamma\delta} + \delta_{\alpha\beta}\varphi \implies \begin{cases} \tau_{xx} = A(x,y)_{,yy} + \phi \\ \tau_{yy} = A(x,y)_{,xx} + \phi \\ \tau_{xy} = -A(x,y)_{,xy} \end{cases} \quad (3)$$

Avec $\phi$ tel que $F_i = -\phi_{,i}$

On peut alors exprimer les contraintes comme le produit de fonctions $S(x,y)$ de x et y, et de valeur nodales $s$ (ou degrés de liberté).

$$\tau = [S(x,y)][s] \quad (4)$$

Le choix de la fonction $A(x,y)$ dépend du type d'élément utilisé. Dans [3], l'élément de Bogner Fox Schmit est utilisé. Il présente quatre degrés de liberté à chaque noeud, qui sont $[A \frac{\partial A}{\partial x} \frac{\partial A}{\partial y} \frac{\partial^2 A}{\partial x \partial y}]^T$. La fonction d'Airy est exprimée par un polynôme d'ordre inférieur à quatre tel que :

$$A(x,y) = \sum_{p,q \le 3} c_{pq} x^p y^q \quad (5)$$

La condition d'équilibre de surface peut être satisfaite en définissant des contacteurs sur chaque face d'un élément et en imposant la continuité des tractions a posteriori par l'utilisation de multiplicateurs de Lagrange [4], ceci augmentant la taille du système à résoudre de manière considérable. Cette condition peut également être satisfaite "a priori". En effet, la continuité du vecteur contrainte entre éléments peut être obtenue en imposant une continuité $C^1$ de la fonction de Airy. La condition d'équilibre de surfaces sera alors vérifiée en imposant l'égalité des degrés de liberté entre éléments.
<figure><img src="image_2.png" /><figcaption>FIG. 1 –</figcaption></figure>

L'élément de Bogner-Fox-Schmit (BFS)

# 3 Procédures numériques

## 3.1 Premier algorithme: algorithme D-PANA

On présente un premier algorithme pour résoudre le problème dans le cas plan. L'idée de cet algorithme, basé sur une technique de relaxation par blocs, consiste en une résolution de deux sous-problèmes : au pas $k$, un problème de frottement à contact donné, et ensuite un problème de contact unilatéral à frottement donné. Nous avons appelé cet algorithme, D-PANA car il est la version duale de l'algorithme classique PANA dû à P. D. Panagiotopoulos [5]. Après condensation [3], on décompose la matrice de compliance $G_h$ qui est symétrique et définie positive en quatre blocs comme suit

$$G_h = \begin{pmatrix} G^{NN} & (G^{TN})^t \\ G^{TN} & G^{TT} \end{pmatrix}, \qquad (6)$$

où $G^{NN}$ (resp. $G^{TT}$) est la sous-matrice correspondant aux degrés de liberté normaux (resp. tangentiels). De façon analogue, le vecteur des déplacements généralisés $g_h$ est décomposé en deux blocs

$$g_h = \begin{pmatrix} g^N \\ g^T \end{pmatrix}, \qquad (7)$$

où $g^N$ (resp. $g^T$) correspond aux degrés de liberté normaux (resp. tangentiels).

Soit $(N^0, T^0)$ donné, l'iteration $k$ de l'algorithme est,

- $T^{k+1}$ est solution du problème de frottement à contact donné suivant

$$\min\{\frac{1}{2}T^tG^{TT}T + T^t(g^T + G^{TN}N^k), |T| \leq -\mu N^k\} \qquad (8)$$

- $N^{k+1}$ est la solution du problème de contact unilatéral contact à frottement donné suivant

$$\min\{\frac{1}{2}N^tG^{NN}N + N^t(g^N + (G^{TN})^tT^{k+1}), N \leq 0\} \qquad (9)$$

Cet algorithme est basé sur une idée simple : si dans l'équation (1), les forces de frottement (resp. de contact) sont données, alors le convexe est fixé et le problème est équivalent à un problème de minimisation sur les forces de contact (resp. de frottement). Notons que les problèmes de minimisation (8)(9) ont des solutions uniques, puisque les fonctionnelles à minimiser sont strictement convexes et que les ensembles dans lesquels on minimise sont convexes et non vides. De plus, On peut résoudre ces problèmes à l'aide de méthodes de minimisation quadratique ou de relaxation efficaces. Il est aussi possible de montrer le résultat suivant :

**Théorème 1** Pour un coefficient de frottement $\mu$ assez petit l'algorithme D-PANA est convergent.

On peut trouver la démonstration de ce résultat dans [6].
## 3.2 Seconde méthode : l’algorithme D-GLT

Cet algorithme est une adaptation de celui proposé dans [7]. L'idée est de minimiser une suite de problèmes en dimension un, c'est à dire, pour la i-ème composante, les autres composantes sont fixées et on a à minimiser un problème convexe à une dimension dans un ensemble convexe

- $N_i \le 0 \quad \forall i \in \{1 \dots n\}$, pour une composante de contact,

- $|T_i| \le -\mu N_i \quad \forall i \in \{1 \dots n\}$, pour une composante de frottement.

L'algorithme est donné par

$(N^0, T^0)$ étant une solution initiale donnée, pour tout $k > 0$ (indice d'itération) et pour $i \in \{1 \dots n\}$:

$$ \alpha = -g_i^N - \sum_{j=1}^{i-1} G_{ij}^{NN} N_j^{k+1} - \sum_{j=i+1}^{n} G_{ij}^{NN} N_j^k - \sum_{j=1}^{i-1} G_{ij}^{TT} T_j^{k+1} - \sum_{j=i}^{n} G_{ij}^{TT} T_j^k \quad (10) $$

$$ N_i^{k+\frac{1}{2}} = \frac{\alpha}{G_{ii}^{NN}} \quad (11) $$

$$ \begin{cases} si \ N_i^{k+\frac{1}{2}} > 0 \ alors \ N_i^{k+1} = 0 \\ si \ N_i^{k+\frac{1}{2}} \le 0 \ alors \ N_i^{k+1} = N_i^{k+\frac{1}{2}} \end{cases} \quad (12) $$

$$ \beta = -g_i^T - \sum_{j=1}^{i} G_{ij}^{NN} N_j^{k+1} - \sum_{j=i+1}^{n} G_{ij}^{NN} N_j^k - \sum_{j=1}^{i-1} G_{ij}^{TT} T_j^{k+1} - \sum_{j=i+1}^{n} G_{ij}^{TT} T_j^k \quad (13) $$

$$ T_i^{k+1} = \frac{\beta}{G_{ii}^{TT}} \quad (14) $$

$$ \begin{cases} si \ T_i^{k+\frac{1}{2}} > -\mu N^{k+1} \ alors \ T_i^{k+1} = -\mu N^{k+1} \\ si \ T_i^{k+\frac{1}{2}} < \mu N^{k+1} \ alors \ T_i^{k+1} = \mu N^{k+1} \\ si \ |T_i^{k+\frac{1}{2}}| \le -\mu N^{k+1} \ alors \ T_i^{k+1} = T_i^{k+\frac{1}{2}} \end{cases} \quad (15) $$

De façon usuelle, cet algorithme est accéléré en introduisant un coefficient de relaxation $\omega$. $N_i^{k+\frac{1}{2}}$ (resp. $T_i^{k+\frac{1}{2}}$) est remplacé par $\omega N_i^{k+\frac{1}{2}} + (1-\omega)N_i^k$ (resp. $\omega T_i^{k+\frac{1}{2}} + (1-\omega)T_i^k$) dans l'algorithme précédent. Notons que, de façon classique, le critère d'arrêt est calculé à partir des variations relatives des forces de contact et de frottement. Ajoutons, pour terminer, que la preuve de convergence de cet algorithme est un problème ouvert.

# 4 Résultats numériques

## 4.1 L'exemple numérique : un "benchmark" classique

Afin de tester et de valider l'algorithme, on a choisi un "benchmark" étudié par le groupe de travail "Validation de codes de calcul" du GRECO "Grandes Déformations et Endommagement" [8]. Le problème est posé en déformations planes. On considère un barreau infini à section rectangulaire, de module d'Young $E = 130$ GPa et de coefficient de Poisson $\nu = 0.2$. Pour des raisons de symétrie, on considère seulement une section carrée. La hauteur est égale à 40mm. La zone de contact est le segment $AB$ sur la figure 2. Le problème a été traité avec différentes valeurs du coefficient de frottement $\mu$ et du chargement ($f$ sur $AD$, $F$ sur $CD$). Les variations des statuts de contact sont données dans la table 1 pour différents jeux de paramètres (les chargements sont en MPa).
<figure><img src="image_3.png" /><figcaption>FIG. 2 – Problème du barreau en contact avec un plan rigide.</figcaption></figure>

<table><caption>TAB. 1 – Comportements de l'interface pour différents paramètres</caption><thead><tr><th>μ</th><th>F</th><th>f</th><th>Gap<br>σ<sub>N</sub> = σ<sub>T</sub> = 0</th><th>Sliding<br>|σ<sub>T</sub>| = μ|σ<sub>N</sub>|</th><th>Stick<br>|σ<sub>T</sub>| < μ|σ<sub>N</sub>|</th></tr></thead><tbody><tr><td>0.2</td><td>-5</td><td>10</td><td>0 mm</td><td>40 mm</td><td>0 mm</td></tr><tr><td>0.2</td><td>-15</td><td>10</td><td>0 mm</td><td>23 mm</td><td>17 mm</td></tr><tr><td>0.2</td><td>-25</td><td>10</td><td>0 mm</td><td>3 mm</td><td>37 mm</td></tr><tr><td>1</td><td>-5</td><td>10</td><td>3 mm</td><td>18 mm</td><td>19 mm</td></tr><tr><td>1</td><td>-5</td><td>15</td><td>3 mm</td><td>25 mm</td><td>12 mm</td></tr></tbody></table>

## 4.2 Concordance avec la formulation primale

Dans ce paragraphe, les résultats numériques sont comparés avec ceux obtenus par la formulation primale [8] afin d'analyser la précision de l'algorithme. Dans [8], le problème est discrétisé par des éléments finis linéaires P1 (triangles à 3 noeuds) et le problème discret est résolu par une méthode de point fixe couplée avec un algorithme de relaxation. Le maillage est non structuré avec 230 noeuds (33 de contact). Dans ce travail, le même problème est traité avec un maillage structuré avec 289 noeuds (17 de contact). Le problème est discrétisé par des éléments de Hermitte bicubiques [4] (quadrangles à 4 noeuds). La figure 3 montre que l'on a une bonne concordance entre les deux formulations en termes de forces de contact. Sur la Figure 4, on montre les déplacements normaux pour les deux formulations, déplacements qui pourraient être un handicap pour l'extension de notre méthode aux grands glissements [3]. Ces figures correspondent au cas $μ = 1$, $F = -5MPa$ and $f = 15MPa$ (5-ième ligne dans la table 1). On note que les algorithmes D-PANA et D-GLT donnent la même solution. Les petites différences entre les courbes peuvent être expliquées de la façon suivante

- les différentes approximations: éléments P1 linéaires ou Hermite bicubiques,

- les différents maillages : triangles avec 230 noeuds ou quadrangles avec 289 noeuds,

- le calcul des déplacements et des forces: par méthode directe ou par dualité.

## Références

[1] M. Jean. Un cadre abstrait pour l'espace vectoriel topologique ordonné $W^{1,p}(\Omega)$ et quelques uns de ses bornés. Technical report, Département de Mathématiques, Université Montpellier 2, 1975. Travaux du séminaire d'analyse convexe, exposé 7.

[2] J.J. Telega. Quasi-static Signorini's contact problem with friction and duality. *International Series of Numerical Mathematics*, 101:199–214, 1991.
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Forces de contact : algorithmes primal et dual.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4 – Déplacements normaux sur la zone de contact : algorithmes primal et dual.</figcaption></figure>

[3] P. Bisegna, F. Lebon, and F. Maceri. Relaxation procedures for solving Signorini-Coulomb contact problems. *Advances in Engineering Software*, 35:595-600, 2004.

[4] V.B. Watwood jr. and B.J. Hartz. An equilibrium stress field model for finite element solutions of two dimensional elastostatic problems. *Int. J. Solids Structures*, 4:857-873, 1968.

[5] P.D. Panagiotopoulos. A nonlinear programming approach to the unilateral contact and friction boundary value problem in the theory of elasticity. *Ingenier Archive*, 44:421-432, 1975.

[6] P. Bisegna, F. Lebon, and F. Maceri. D-PANA: a convergent block-relaxation solution method for the discretized dual formulation of the Signorini-Coulomb contact problem. *Comptes Rendus Académie des Sciences, Paris, Série I*, 333:1053-1058, 2001.

[7] R. Glowinski, J. L. Lions, and R. Trémolières. *Analyse numérique des inéquations variationnelles*. Dunod, 1975.

[8] M. Raous, P. Chabrand, and F. Lebon. Numerical methods for solving unilateral contact problem with friction. *Journal of Theoretical and Applied Mechanics*, 7:111-128, 1988.