# Extension de la méthode du bi-potentiel à l'analyse des problèmes d'impact

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

► To cite this version:

Benoit Magnain, Zhi-Qiang Feng, Jean-Michel Cros. Extension de la méthode du bi-potentiel à l'analyse des problèmes d'impact. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813055

HAL Id: hal-01813055

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

Submitted on 12 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 CC0 1.0 - Universal - International License
# Extension de la méthode du bi-potentiel à l'analyse des problèmes d'impact

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

Laboratoire de Mécanique et d'Énergétique d'Évry
Université d'Évry - Val d'Essonne
40, rue du Pelvoux
F-91020 Évry Cedex
{magnain,feng,cros}@iup.univ-evry.fr

ABSTRACT. The bi-potential method has been successfully applied for the modelling of frictional contact problems in static cases. This paper presents the extension of this method for dynamic analysis of impact problems. A first order algorithm is applied for the numerical integration of the time-discretized equation of motion. The solution algorithm is simple and efficient. In order to validate the proposed method, numerical results for a longitudinal impact between two elastic bars are presented.

MOTS-CLÉS : Contact, Frottement, Dynamique, Méthode du bi-potentiel.

KEYWORDS: Contact, Friction, Dynamic, Bi-potentiel method.

<footnote>RÉSUMÉ. La méthode du bi-potentiel a été appliquée avec succès à la modélisation des problèmes de contact avec frottement pour les cas statiques. Cet article présente l'extension de cette méthode aux problèmes d'impact en dynamique. Un algorithme du premier ordre est adopté pour l'intégration des équations du mouvement. L'algorithme ainsi obtenu est simple et efficace. Afin de valider la méthode proposée, les résultats de la simulation numérique de l'impact longitudinal entre deux barres élastiques sont présentés.</footnote>
# 1. Introduction

Les problèmes de contact avec frottement interviennent dans un nombre considérable d'applications industrielles. Dans le domaine des méthodes numériques, de nombreux travaux ont été menés afin de développer de nouveaux algorithmes de calcul capables de traiter efficacement ce type de problème. Actuellement, les codes de calcul industriels (ANSYS, PAM-CRASH...) sont presque exclusivement basés sur les méthodes de pénalisation et de régularisation. Pourtant, pour ce type de méthodes, les conditions de contact et les lois de frottement ne sont pas satisfaites exactement. De plus, les paramètres de régularisation sont parfois délicats à déterminer. Pour éviter ces difficultés, une méthode du Lagrangien augmenté a été développée par Alart et Curnier [ALA 91]. Cette méthode repose sur l'utilisation d'un algorithme de Newton généralisé pour déterminer simultanément toutes les inconnues du problème (déplacement et réactions). De même, De Saxcé et Feng [DES 91] ont proposé une méthode du bi-potentiel dans laquelle le problème du 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. Cette dernière méthode n'utilise qu'un seul principe variationnel sur le déplacement et une seule inégalité sur le contact unilatéral et le frottement. Cette formulation du Lagrangien augmenté a été appliquée avec succès aux problèmes de contact avec frottement pour des cas quasi statiques. Dans ce travail, nous proposons une extension de cette méthode à l'analyse des problèmes d'impact entre plusieurs corps déformables. Le plus souvent, la discrétisation temporelle des équations du mouvement se fait à l'aide de schémas d'intégration du second ordre tels que les schémas de Newmark et Wilson. Ces schémas supposent que l'accélération est constante ou qu'elle varie de manière linéaire à l'intérieur d'un pas de temps. Or, dans les problèmes d'impact, la vitesse et l'accélération ne sont pas continues au moment du choc en raison des changements brutaux des conditions de contact. C'est pourquoi il apparaît judicieux d'adopter un schéma du premier ordre comme celui développé par Jean [JEA 99], dans lequel le calcul de l'accélération n'intervient pas. Le modèle ainsi développé combine la méthode du bi-potentiel pour la résolution du problème de contact et un schéma du premier ordre pour la discrétisation temporelle. Afin de valider la méthode développée, nous proposons de l'appliquer à l'étude d'un impact longitudinal entre deux barres élastiques. Ce cas simple, déjà traité par Hu [HU 97], permet en effet de comparer les résultats numériques obtenus avec la solution analytique fournie par la théorie 1D. De plus, il met en évidence une propriété intéressante de la méthode développée : le principe de conservation d'énergie est satisfait dans le cas d'un contact sans frottement et cela sans avoir recours à une technique de régularisation. D'autres exemples plus complexes sont présentés dans un papier compagnon [FEN 05]. En particulier, des cas d'impact avec frottement entre plusieurs corps déformables y sont traités avec succès.
## 2. Méthode du bi-potentiel

### Contact unilatéral et modèle de frottement de Coulomb

Soit $\mathbf{r}$ la réaction de contact et $\dot{\mathbf{u}}$ la vitesse relative des deux corps en contact. De manière classique, le contact unilatéral est décrit par les conditions de Signorini définies sur l'ensemble des points en contact. Dans le cadre de la dynamique, ces conditions s'écrivent dans le repère local du contact en terme de vitesse relative normale $\dot{u}_n$ et de réaction normale $r_n$ :

$$ \dot{u}_n \ge 0, \quad r_n \ge 0 \text{ et } r_n \dot{u}_n = 0 \qquad (1) $$

Dans le cas où l'on adopte le modèle de frottement de Coulomb, l'ensemble des réactions admissibles est représenté par le convexe $K_\mu$ défini par :

$$ K_{\mu} = \{ \mathbf{r} \in \mathbb{R}^3 \text{ tel que } |\mathbf{r}_{\mathbf{t}}| \le \mu r_n \} \qquad (2) $$

De Saixé et Feng [DES 98] ont montré que la loi de contact peut s'écrire :

$$ -(\dot{\mathbf{u}}_{\mathbf{t}} + (\dot{u}_n + \mu|\dot{\mathbf{u}}_{\mathbf{t}}|)\mathbf{n}) \in \partial \bigcup_{K_{\mu}} (\mathbf{r}) \quad \text{où} \quad \bigcup_{K_{\mu}} (\mathbf{r}) = \begin{cases} 0 & \text{si } \mathbf{r} \in K_{\mu} \\ +\infty & \text{sinon} \end{cases} \qquad (3) $$

Ils définissent alors le bi-potentiel de contact suivant :

$$ 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 \text{avec } \mathbb{R}_- = ]-\infty, 0] \qquad (4) $$

### Algorithme local

Dans le but d'éviter les potentiels non différentiables qui apparaissent en mécanique non linéaire, comme c'est le cas en mécanique du contact, il est souhaitable d'utiliser la méthode du Lagrangien augmenté. Pour le bi-potentiel de contact défini par (4), et sous la condition que $\dot{u}_n \ge 0$ et $\mathbf{r} \in K_\mu$, nous définissons les réactions de contact augmentées $\boldsymbol{\tau}$ par :

$$ \boldsymbol{\tau} = \mathbf{r} - \rho(\dot{\mathbf{u}}_t + (\dot{u}_n + \mu|\dot{\mathbf{u}}_t|)\mathbf{n}) \qquad (5) $$

avec $\rho > 0$ déterminé de manière à assurer la convergence du calcul et l'on montre que $\mathbf{r}$ est la projection de $\boldsymbol{\tau}$ sur le cône de Coulomb :

$$ \mathbf{r} = \operatorname{proj}(\boldsymbol{\tau}, K_{\mu}) \qquad (6) $$

La résolution numérique de l'équation (6) se fait à l'aide de l'algorithme d'Uzawa. Ce dernier conduit à un processus itératif incluant un schéma de type prédiction-correction :

$$ \left\{ \begin{array}{l} \text{prédiction } \boldsymbol{\tau}^{i+1} = \mathbf{r}^i - \rho^i (\dot{\mathbf{u}}_t^i + (\dot{u}_n^i + \mu|\dot{\mathbf{u}}_t^i|)\mathbf{n}) \\ \text{correction } \mathbf{r}^{i+1} = \operatorname{proj}(\boldsymbol{\tau}^{i+1}, K_\mu) \end{array} \right. \qquad (7) $$

Un point essentiel de la méthode du bi-potentiel est que la correction peut être évaluée analytiquement :

$$ \left\{ \begin{array}{ll} \text{si} & \mu|\tau_t^{i+1}| < -\tau_n^{i+1} \text{ alors } \mathbf{r}^{i+1} = 0 \qquad (\text{pas de contact}) \\ \text{sinon si} & |\tau_t^{i+1}| < \mu \tau_n^{i+1} \text{ alors } \mathbf{r}^{i+1} = \tau^{i+1} \qquad (\text{adhérence}) \\ \text{sinon} & \mathbf{r}^{i+1} = \tau^{i+1} - \frac{(|\tau^{i+1}| - \mu \tau_n^{i+1})}{(1+\mu^2)} \left( \frac{\tau^{i+1}}{|\tau^{i+1}|} + \mu \mathbf{n} \right) \qquad (\text{glissement}) \end{array} \right. \qquad (8) $$
### 3. Algorithme du premier ordre

Après discrétisation spatiale par la méthode des éléments finis, un problème de contact en dynamique et dans le cas non-linéaire s'écrit sous la forme de l'équation d'équilibre suivante :

$$ M \ddot{u} = F + R_c, \quad \text{où } F = F_{ext} + F_{int} - C \dot{u} \quad (9) $$

Les vecteurs $F_{ext}$, $F_{int}$ et $R_c$ représentent respectivement les forces externes, internes, et de contact. M est la matrice de masse du système et C la matrice d'amortissement. L'intégration de l'équation (9) est réalisée par un schéma du premier ordre [JEA 99]. Les équations ne font alors plus intervenir l'accélération qui pose des difficultés au moment de l'impact. La linéarisation du vecteur force donne :

$$ F^{i+1} = F_{int}^i + \frac{\partial F}{\partial u}(u^{i+1} - u^i) + \frac{\partial F}{\partial \dot{u}}(\dot{u}^{i+1} - \dot{u}^i) = F_{int}^i - K^i \Delta u - C^i \Delta \dot{u} \quad (10) $$

Finalement, nous obtenons la forme récursive en déplacement suivante ($0 \le \xi \le 1$ et $0 \le \theta \le 1$) :

$$ \left\{ \begin{array}{l} \bar{K}^i \Delta u = \bar{F}^i + \bar{F}_{acc}^i + \bar{R}_c^{i+1} \\ \mathbf{u}^{i+1} = \mathbf{u}^i + \Delta \mathbf{u} \end{array} \right. \quad (11) $$

où $\bar{K}^i$, $\bar{F}_{acc}^i$ et $\bar{F}^i$ sont calculés par :

$$ \left\{ \begin{array}{l} \bar{K}^i = \xi K^i + \frac{\xi}{\theta \Delta t} C^i + \frac{1}{\theta \Delta t^2} M^i \\ \bar{F}_{acc}^i = -\frac{1}{\theta \Delta t^2} M^i \{\mathbf{u}^i - \mathbf{u}^t - \Delta t \dot{\mathbf{u}}^t\} \\ \bar{F}^i = (1-\xi)(\bar{F}_{int}^t + \bar{F}_{ext}^t) + \xi (\bar{F}_{int}^i + \bar{F}_{ext}^{t+\Delta t}) \end{array} \right. \quad (12) $$

À 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} \left(\mathbf{u}^{t+\Delta t} - \mathbf{u}^t\right) \quad (13) $$

Pour $\theta = 0, 5$, 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.

### 4. Résultats numériques

L'algorithme présenté précédemment a été implanté et testé dans le code éléments finis FER/Impact. Afin d'illustrer l'efficacité et la robustesse de la méthode, nous considérons un exemple d'application appartenant à la catégorie des problèmes homogènes de Neumann (aucun déplacement et aucune force extérieure ne sont impossés). Il s'agit du problème, déjà largement étudié, de l'impact longitudinal entre deux barres élastiques en 2D. La configuration géométrique du problème est reportée sur la figure 1. Pour ce problème adimensionné, les valeurs numériques utilisées sont les suivantes : $L = 10, S = 1, E = 1000, \nu = 0$ et $\rho = 0,001$. La distance initiale entre les
deux barres est égale à 0,02. Chacune des deux barres est maillée par l'intermédiaire de 20 éléments à 4 noeuds. Les vitesses initiales des deux barres sont respectivement :
$$V_0^{(1)} = -V_0^{(2)} = 1.$$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Impact longitudinal entre deux barres élastiques</figcaption></figure>

Les valeurs choisies pour les paramètres de l'algorithm du premier ordre sont les suivantes : $\theta = \xi = 0,5$ et $\Delta t = 10^{-5}$. De plus, nous considérons que l'amortissement du matériau est négligeable dans ce problème i.e. $C = 0$ dans Eqs.(9-12). Pour chaque résultat présenté dans la suite, le résultat analytique obtenu à l'aide de la théorie 1D est proposé afin de vérifier la précision des calculs. L'évolution du déplacement et de la vitesse du point C est tracée sur les figures 2 et 3. Ces deux graphiques montrent la qualité des résultats obtenus à l'aide de la méthode proposée. Les oscillations de la solution numérique (figure 3) traduisent les vibrations de la barre suite à l'impact phénomène que ne peut représenter le modèle analytique. Ces oscillations indiquent la présence d'une énergie de déformation élastique (même faible) dans la barre. La figure 4 propose la comparaison entre la force de réaction calculée à l'aide de l'algorithm et celle proposée par la théorie. Là encore, les courbes montrent que la méthode proposée permet d'obtenir une grande précision sur les résultats. L'évolution de l'énergie totale ($E_t$) du système (barre(1)+barre(2)) est tracée sur la figure 5 ainsi que les évolutions des énergies cinétiques et élastique du système ($E_c$ et $E_e$). On peut constater que l'énergie totale du système est effectivement conservée.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Déplacement du point C</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Vitesse du point C</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. Force de réaction au point C</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Évolution de l'énergie</figcaption></figure>

**5. Conclusion**

Dans cet article, nous avons présenté les récents développements apportés à la méthode du bi-potentiel destinés à l'étude des problèmes d'impact. Une description détaillée de l'algorithme obtenu est proposée. Les résultats numériques montrent que la méthode est robuste et efficace. Par ailleurs, ils permettent de penser que l'algorithme préserve l'énergie totale du système pour les cas de contact sans frottement.

**6. 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.

[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 05] FENG Z.-Q., MAGNAIN B., CROS J.-M., « FER/Impact : Logiciel de simulation numérique des problèmes d'impact », *7ème Colloque National en Calcul des Structures*, 2005.

[HU 97] HU N., « A solution method for dynamic contact problems », *Computers and Structures*, vol. 63, n° 6, 1997, p. 1053-1063.

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

[SIM 91] SIMO J., WONG 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.