# Résolution du problème de contact mécanique frottant : méthode du lagrangien augmenté adapté

Philippe Bussetta, Daniel Marceau, Jean-Philippe Ponthot

► To cite this version:

Philippe Bussetta, Daniel Marceau, Jean-Philippe Ponthot. Résolution du problème de contact mécanique frottant : méthode du lagrangien augmenté adapté. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422251

HAL Id: hal-01422251

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

Submitted on 24 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
# Résolution du problème de contact mécanique frottant :
méthode du lagrangien augmenté adapté

P. Bussetta¹, D. Marceau¹, J.-P. Ponthot²

¹ Université du Québec à Chicoutimi (Canada)
555, boulevard de l'Université Chicoutimi (Québec, Canada) G7H 2B1
{Philippe_Bussetta,Daniel_Marceau}@uqac.ca

² Université de Liège (Belgique)
LTAS MN²L Institut de mécanique et de génie civil, Bâtiment B52/3
Chemin des Chevreuils, 1 B-4000 Liege Belgique
jp.ponthot@ulg.ac.be

**Résumé** — Cet article présente une nouvelle méthode de résolution du problème de contact mécanique frottant, celle du lagrangien augmenté adapté. Contrairement aux méthodes les plus utilisées (pénalisation et lagrangien augmenté), l'utilisateur ne doit pas déterminer les valeurs des coefficients de pénalisation. Cette méthode est performante aussi bien avec des matériaux élastiques qu'élasto-plastiques. Le problème de Hertz est présenté. Cette méthode est aussi utilisée dans le contexte multi-physique.

**Mots clés** — Contact mécanique, frottement, Lagrangien augmenté, Lagrangien augmenté adapté, pénalisation, pénalité adaptative.

## 1 Introduction

Le contact mécanique demeure, encore aujourd'hui, le problème de mécanique des solides qui présente les non-linéarités les plus difficiles à prendre en compte. Les méthodes les plus utilisées pour formuler un système d'équations sont celles de pénalisation et du lagrangien augmenté. Ces méthodes demeurent difficiles à utiliser en raison de la difficulté d'identification des valeurs des coefficients de pénalisation. Afin de pallier les carences de ces méthodes, une nouvelle approche est proposée, celle du « lagrangien augmenté adapté ». Cette méthode est basée sur celle du lagrangien augmenté jumelée à une adaptation de la pénalité inspirée de celle proposée par Chamoret [2]. Avec cette nouvelle méthode, l'utilisateur n'est plus tenu de choisir les valeurs des coefficients de pénalisation. La fiabilité et la rapidité de cette approche sont prouvées avec la résolution du problème de Hertz. Ces travaux font partie d'un projet de recherche portant sur la modélisation et la résolution du problème de contact mécanique dans un contexte multi-physique [1].

## 2 Équations du problème

Le problème de contact mécanique entre deux solides déformables est traité. La partie de leur frontière pouvant rentrer en contact est notée dans la configuration d'origine $\Gamma_c^1$ pour le solide 1 et $\Gamma_c^2$ pour le solide 2. Le principe des travaux virtuels peut s'écrire sous la forme suivante :

$$ \mathcal{P}^{\text{int}}(u, \delta u) + \mathcal{P}^{\text{ext}}(u, \delta u) + \mathcal{P}^{\text{c}}(u, \delta u) = 0 \quad \text{avec} \quad \mathcal{P}^{\text{c}}(u, \delta u) = \int_{\Gamma_c^1} (\delta\xi^\beta t_{t\beta} + \delta g t_n) d\Gamma_c^1 $$
u désigne le champ de déplacement réel et $\delta u$ le virtuel. $\mathcal{P}^{\text{int}}(u, \delta u)$ représente le travail virtuel produit par les forces internes et $\mathcal{P}^{\text{ext}}(u, \delta u)$ celui des forces externes à l'exception des forces de contact. Ces dernières sont incluses dans le terme $\mathcal{P}^c(u, \delta u)$ ($\beta = 1$ ou $2$). $\delta\xi^\beta$ représente le champ de déplacement virtuel tangentiel et $\delta g$ celui normal sur $\Gamma_c^1$ (pour plus d'informations voir [3]). La contrainte de contact est décomposée en une partie normale ($t_n$) et une tangentielle ($t_{t_\beta}$) [1]. Les utilisées pour modéliser le contact frottant sont, celle du contact unilatéral (relations (1)) et celle du frottement de Coulomb (relations (2)).

$$g \le 0 \qquad t_n \ge 0 \qquad gt_n = 0 \qquad (1)$$

$$\Phi = \|t_t\| - \mu t_n \le 0 \qquad \zeta \ge 0 \qquad \Phi \zeta = 0 \qquad (2)$$

$\zeta$ est le taux de glissement absolu, $g$ représente la distance entre les deux solides et $\mu$ est le coefficient de frottement de Coulomb.

# 3 Méthodes de résolution

## 3.1 Méthode de pénalisation

L'introduction de coefficient de pénalisation normale ($\varepsilon_n$) et tangentielle ($\varepsilon_t$) permet d'obtenir des relations continues liant les contraintes de contact aux déplacements :

$$\begin{aligned} t_n &= \varepsilon_n \langle g \rangle & \Phi &= \|t_t\| - \mu t_n \le 0 \\ v_t - \zeta \frac{t_t}{\|t_t\|} &= \frac{1}{\varepsilon_t} \mathcal{L}_v t_t \quad \text{avec} \quad \zeta \ge 0 & \Phi \zeta &= 0 \end{aligned}$$

L'opérateur $\mathcal{L}_v$ désigne la dérivée de Lie et $\langle \rangle$ la partie positive de son opérande. Le champ de déplacement ($u$) est obtenu en suivant l'algorithme de Newton-Raphson.

Le choix des coefficients de pénalisation est primordial. Des valeurs trop faibles de ces coefficients engendrent des erreurs sur les conditions de contact qui n'ont aucun sens physique. Alors que des valeurs trop grandes provoquent des oscillations numériques pouvant empêcher la convergence. Une adaptation de la pénalité a été proposée par Chamoret [2], mais uniquement au niveau de la pénalisation normale et elle n'est pas bien adaptée aux problèmes élasto-plastiques [1].

## 3.2 Méthode du lagrangien augmenté

La résolution est effectuée de la même manière qu'avec la méthode de pénalisation, mais les contraintes de contact sont définies par les relations suivantes :

$$\begin{aligned} t_n &= \langle \lambda_n + \varepsilon_n g \rangle & \Phi &= \|t_t\| - \mu t_n \le 0 \\ v_t - \zeta \frac{t_t}{\|t_t\|} &= \frac{1}{\varepsilon_t} (\mathcal{L}_v t_t - \mathcal{L}_v \lambda_t) \quad \text{avec} \quad \zeta \ge 0 & \Phi \zeta &= 0 \end{aligned}$$

Les multiplicateurs de Lagrange ($\lambda$) sont décomposés en une partie normale ($\lambda_n$) et une tangentielle ($\lambda_t$). Ils sont initialisés avec les valeurs des contraintes de contact du pas de temps précédent. Ils sont mis à jour suivant les relations (3) lors de la convergence de l'algorithme de Newton-Raphson, si l'interpénétration ou le glissement réversible sont supérieurs aux limites.

$$\lambda_n^{k+1} = t_n^{k+1} = \langle \lambda_n^k + \varepsilon_n g \rangle \qquad \lambda_t^{k+1} = t_t^{k+1} \qquad (3)$$

Les augmentations des multiplicateurs de Lagrange autorisent l'utilisation de valeurs plus faibles des coefficients de pénalisation. Néanmoins, une valeur trop faible de ces coefficients peut engendrer des temps de calcul très importants en raison de la faible vitesse de convergence.
### 3.3 Méthode du lagrangien augmenté adapté

La méthode du lagrangien augmenté adapté est basée sur celle du lagrangien augmenté et sur une adaptation de la pénalisation. La résolution est effectuée de la même manière qu'avec la méthode du lagrangien augmenté. Mais, la valeur des coefficients de pénalisation et des contraintes de contact sont calculées différemment.

<caption>**Algorithm 1** Calcul de la valeur de la contrainte normale de contact</caption>

**Require:** $\varepsilon_n$, $g_i$ et $g_{i-1}$.
**if** $g_i \times g_{i-1} < 0$ **then**
    **if** deuxième changement consécutif **then**
        $Evalt_n = 0$
    **else**
        $Evalt_n = \varepsilon_n g_{i-1}$
    **end if**
**end if**
**if** $g_i < 0$ ou $t_n(i-1) > 0$ **then**
    $\varepsilon_n$ = adaptation de $\varepsilon_n$
    **else**
        $\varepsilon_n$ = valeur initiale de $\varepsilon_n$
    **end if**
    $t_n = < Evalt_n + \varepsilon_n g_i + \lambda_n >$

La contrainte normale de contact est calculée de la même manière qu'avec la méthode du lagrangien augmenté. Mais un terme est ajouté afin de mieux l'évaluer lorsque la valeur de l'interpénétration change de signe (voir l'algorithme 1). Cette procédure de calcul est basée sur le principe suivant : si la valeur de l'interpénétration change de signe ($g_i \times g_{i-1} < 0$), cela signifie que la correction apportée est trop importante. Cette dernière ($g_{i-1}\varepsilon_n$) est ajoutée à la composante normale du multiplicateur de Lagrange afin de réduire les oscillations numériques. Par contre, lorsque ce phénomène se répète plusieurs fois consécutives, aucun historique des corrections apportées par la pénalisation n'est conservé car la valeur trop importante du coefficient de pénalisation normale provoque ces instabilités (voir l'algorithme 1).

Lors de l'entrée en contact, une valeur initiale du coefficient de pénalisation normale est impossée. Par la suite, cette valeur est ajustée en fonction de l'historique des conditions de contact (voir l'algorithme 2). Trois cas de figures sont envisageables, soit la valeur de l'interpénétration change de signe ($g_i \times g_{i-1} < 0$), soit la valeur absolue de l'interpénétration est supérieure à la limite maximale ($|g_i| > g_{max}$), soit elle est inférieure à la limite maximale ($|g_i| < g_{max}$). Dans le dernier cas, la valeur du coefficient de pénalisation normale est identique.

Dans le cas où $g_i \times g_{i-1}$ est négatif, la valeur de $\varepsilon_n$ est ajustée afin de rendre la valeur de l'interpénétration la plus petite possible en valeur absolue sans qu'elle change à nouveau de signe.

Dans le cas où $|g_i|$ est supérieur à $g_{max}$, l'adaptation de $\varepsilon_n$ doit permettre d'augmenter la contribution de la pénalisation dans la valeur de la contrainte normale sans que l'interpénétration change de signe. Afin d'optimiser le temps de calcul plusieurs sous-cas sont envisagés. Si la valeur de l'interpénétration varie trop, alors afin d'éviter les oscillations numériques, le coefficient de pénalisation normale est multiplié par deux. Si la valeur de l'interpénétration est quasiment constante, alors afin de réduire le temps de calcul, le coefficient de pénalisation normale est adapté deux fois (le coefficient multiplicatif est élevé au carré) dans la limite où cette augmentation n'est pas trop importante. Si la valeur absolue de l'interpénétration augmente, alors afin de la ramener le plus rapidement possible sous la limite sans engendrer d'oscillations numériques, l'adaptation de $\varepsilon_n$ permet de doubler la contribution de la pénalisation dans la valeur de la contrainte normale. Si non, la fonction par laquelle le coefficient de pénalisation normale est multiplié, a une croissance
plus lente que celle utilisée avec la méthode d'adaptation de la pénalité (linéaire [2]), ce qui réduit les risques d'oscillations numériques.

<caption>**Algorithm 2** Adaptation du coefficient de pénalisation normale</caption>

**Require:** $\varepsilon_n, g_i$ et $g_{i-1}$.
**if** $g_i \times g_{i-1} < 0$ **then**
    **if** $g_{i-1} > g_{max}$ **then**
        $\varepsilon_n = |(\varepsilon_n g_{i-1})/g_i \times (|g_i| + g_{max})/(g_i - g_{i-1})|$
    **else**
        $\varepsilon_n = |\varepsilon_n g_{i-1}|/(10g_i)$
    **end if**
    **else if** $g_i > g_{max}$ **then**
        **if** $|g_i - g_{i-1}| > max(g_i/10; g_{i-1}/10; 5g_{max})$ **then**
            $\varepsilon_n = 2\varepsilon_n$
        **else if** $|0.99g_{i-1}| < |g_i| < |1.01g_{i-1}| < 10g_{max}$ **then**
            $\varepsilon_n = \varepsilon_n(\sqrt{(|g_i|/g_{max} - 1) + 1})^2$
        **else if** $|g_i| > |1.01g_{i-1}|$ **then**
            $\varepsilon_n = 2\varepsilon_n(g_{i-1}/g_i)$
        **else**
            $\varepsilon_n = \varepsilon_n(\sqrt{(|g_i|/g_{max} - 1) + 1})$
    **end if**
    **else**
        $\varepsilon_n = \varepsilon_n$
    **end if**

Lors de l'augmentation des multiplicateurs de Lagrange, la valeur maximale de l'interpénétration pour l'adaptation de la pénalité est remplacée par la valeur de l'interpénétration maximale divisée par dix. Ainsi, l'utilisateur peut imposer une interpénétration maximale aussi petite qu'il le désire. Pour cela, il doit soit choisir deux de ces trois valeurs : un nombre de cycle d'augmentation des multiplicateurs de Lagrange, la valeur maximale de l'interpénétration pour la solution ou celle lors de la première augmentation des multiplicateurs de Lagrange.

<caption>**Algorithm 3** Calcul de la valeur du coefficient de pénalisation tangentielle</caption>

**Require:** $\varepsilon_t$, GLRe = Valeur absolue du glissement réversible.
**if** (pas de contact à l'itération $i-1$) ou (Première itération après l'augmentation des multiplicateurs de Lagrange et GLRe < Tol) **then**
    $\varepsilon_t = \mu t_n / (10 \text{ GLRe})$
**else if** retournement du sens de glissement **then**
    $\varepsilon_t = \varepsilon_t / 2$
**end if**

La valeur du coefficient de pénalisation est basée sur une étude empirique qui montre qu'une surévaluation de la contrainte tangentielle perturbe la résolution. Lors de l'entrée en contact, la valeur du coefficient de pénalisation tangentielle permet d'obtenir une contrainte tangentielle égale à dix pour cent de celle de glissement. La valeur du coefficient de pénalisation tangentielle est ajustée pendant l'augmentation des multiplicateurs de Lagrange (voir l'algorithme 3). De plus, la valeur du coefficient de pénalisation tangentielle est divisée par deux si un retournement du sens de glissement se produit entre deux itérations.
# 4 Résultats numériques

Le problème de Hertz entre deux cylindres identiques est traité en deux dimensions [4]. Dans le cas où le matériau est élastique, l'hypothèse des déformations planes est retenue; alors que, dans le cas élasto-plastique, celle des contraintes planes est utilisée [1]. Deux quart de cylindre sont modélisés, ils sont maillés de manière identique. Trois maillages différents sont utilisés (noté maillages n où n est le nombre d'éléments dans la zone de contact, n=2, 4 ou 6).

<table><caption>Tableau 1 – Nombre d'itérations (calcul fait en deux pas de temps)</caption><thead><tr><th rowspan="2"></th><th rowspan="2">ε<sub>n</sub><sup>0</sup></th><th colspan="3">élastique</th><th colspan="3">élasto-plastique</th></tr><tr><th>2</th><th>4</th><th>6</th><th>2</th><th>4</th><th>6</th></tr></thead><tbody><tr><td>pénalisation</td><td>10<sup>6</sup></td><td>4-3</td><td>4-3</td><td>6-4</td><td>5-7</td><td>5-6</td><td>os ct</td></tr><tr><td>lagrangien</td><td>2</td><td>14-14</td><td>15-15</td><td>15-15</td><td>15-14</td><td>15-15</td><td>15-17</td></tr><tr><td>augmenté</td><td>2 × 10<sup>2</sup></td><td>14-12</td><td>12-12</td><td>14-15</td><td>14-15</td><td>12-13</td><td>13-12</td></tr><tr><td>adapté</td><td>2 × 10<sup>4</sup></td><td>7-7</td><td>8-8</td><td>9-8</td><td>8-9</td><td>7-9</td><td>9-7</td></tr></tbody></table>

$$ε_n^0 = \text{valeur initiale de } ε_n$$

Le tableau 1 et les figures 1 et 2 permettent de constater qu'avec la méthode du lagrangien augmenté adapté la valeur initiale du coefficient de pénalisation normale n'a pas d'influence sur la qualité de la solution et que très peu sur le temps de calcule. La méthode du lagrangien augmenté adapté nécessite plus de temps de calcul que celle de pénalisation. Mais, cette dernière ne permet pas de résoudre le problème élasto-plastique pour le maillage 6 à cause d'oscillations numériques des statuts de contact. De plus, avec la méthode de pénalisation la valeur des coefficients de pénalisation a une grande influence sur la valeur des contraintes de contact (figures 1(a) et 1(b)).

<figure><img src="image_2.png" /><figcaption>Figure 1 – Contrainte normale de contact (pénalisation, maillage 4)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2 – Contrainte normale de contact (lagrangien augmenté adapté, maillage 4)</figcaption></figure>

**5 Conclusion**

Cet article présente une nouvelle méthode de résolution du problème de contact mécanique frottant, la méthode du lagrangien augmenté adapté. Cette méthode possède les avantages des méthodes les plus utilisées. Elle permet de contrôler l'interpénétration et présente la stabilité de la méthode du lagrangien augmenté et la vitesse de celle de la pénalisation. De plus, l'utilisateur ne doit pas déterminer la valeur numérique des coefficients de pénalisation. Cette méthode permet d'obtenir plus facilement, une meilleure évaluation des contraintes de contact. Elle est aussi utilisée dans le contexte multi-physique pour la résolution de problème industriel thermoélectro-mécanique [1]. La méthode du lagrangien augmenté adapté peut facilement être utilisée pour prendre en compte les déformations des aspérités des surfaces de contact.

**Références**

[1] P. Bussetta. *Modélisation et résolution du problème de contact mécanique et son application dans un contexte multi-physique*. Thèse de doctorat, Université du Québec à Chicoutimi, Canada 2009

[2] D. Chamoret *Modélisation du contact : nouvelles approches numériques* Thèse de doctorat, École Centrale de Lyon, France 2002

[3] D. Marceau. *Modélisation du contact tridimensionnel avec frottement en grande transformation et son application à l'étude des dispositifs d'ancrage multitorons* Thèse de doctorat, Université Laval (Québec), Canada 2001

[4] S. Timoshenko, J.N. Goodier. *Theory of Elasticity*. Second Edition, Mc Graw-Hill Book Company, Inc, New York 1951.