# Polarisation des forces de friction dans un système multi-contacts

Serge Dumont, Jérôme Fortin, Youssef Ouafik

► To cite this version:

Serge Dumont, Jérôme Fortin, Youssef Ouafik. Polarisation des forces de friction dans un système multi-contacts. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413708

HAL Id: hal-01413708

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

Submitted on 10 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
# Polarisation des forces de friction dans un système multi-contacts

S. Dumont¹, J. Fortin², Y. Ouafik³

¹ LAMFA UPJV - CNRS UMR 6140
33 rue Saint-Leu, 80 000 Amiens, France
serge.dumont@u-picardie.fr

² INSSET UPJV - LTI EA 3899
48 rue Raspail, 02100 Saint-Quentin, France
jerome.fortin@u-picardie.fr

³ CEMEF - CNRS UMR 7635
BP 207, 06 904 Sofia Antipolis cedex, France
youssef.ouafik@mines-paristech.fr

**Résumé** — Le contexte du travail proposé ici consiste à essayer d'identifier des lois de comportement homogénéisées rendant compte du comportement global des systèmes multi-contacts. A l'heure actuelle, ces systèmes sont généralement analysés avec la Mécanique des Milieux Continus et la micro-mécanique ou approche multi-échelle. Cette approche se distingue de l'approche phénoménologique classiquement utilisée pour la modélisation du comportement des matériaux solides qui est basée sur des formulations mathématiques développées dans le cadre de la thermodynamique et dont les constantes sont définies à partir de résultats d'essais de laboratoire.

**Mots clés** — Méthode des Eléments Discrets ; Méthode multi-échelle ; Frottement.

## 1 Introduction

Le contexte du travail proposé ici consiste à essayer d'identifier des lois de comportement homogénéisées rendant compte du comportement global des systèmes multi-contacts. A l'heure actuelle, ces systèmes sont généralement analysés avec la Mécanique des Milieux Continus et la micro-mécanique ou approche multi-échelle. Cette approche se distingue de l'approche phénoménologique classiquement utilisée pour la modélisation du comportement des matériaux solides qui est basée sur des formulations mathématiques développées dans le cadre de la thermodynamique et dont les constantes sont définies à partir de résultats d'essais de laboratoire. Le manque de physique de base dans ces formulations mathématiques génère des modèles souvent complexes et difficiles à identifier. L'approche multi-échelle paraît bien adaptée pour palier ces difficultés. Cette étude a pour but de quantifier à l'aide de la méthode des Eléments Discrets [6, 7] les phénomènes de polarisation des efforts de contact.
## 2 Approche par Eléments Discrets

Dans les domaines où la collection de corps entre lesquels des liaisons unilatérales, usuellement affectées du frottement de Coulomb, sont susceptibles de s'établir ou de se rompre, ces équations de contraintes conduisent à un problème de complémentarité non linéaire qui ne peut être résolu par une méthode de programmation linéaire.

A chaque paire de particules $\Omega_i$ et $\Omega_j$ candidates au contact, est associé un référentiel local dont les axes sont orientés suivant les deux vecteurs unitaires $\mathbf{n}$ et $\mathbf{t}$, respectivement normal et tangent au plan de contact. La normale $\mathbf{n}$ est dirigée de $\Omega_j$ vers $\Omega_i$. Les variables mises en dualité sont $\dot{\mathbf{u}}$, la vitesse relative locale de $\Omega_i$ par rapport à $\Omega_j$, et la réaction de contact $\mathbf{r}$ de $\Omega_j$ sur $\Omega_i$. Elles se décomposent dans la base locale par :

$$ \dot{\mathbf{u}} = \dot{\mathbf{u}}_t + \dot{u}_n \cdot \mathbf{n}, \quad \mathbf{r} = \mathbf{r}_t + r_n \cdot \mathbf{n}, \qquad (1) $$

où $\dot{u}_n$ est la vitesse de séparation normale, $\dot{\mathbf{u}}_t$ la vitesse de glissement, $r_n$ la pression de contact et $\mathbf{r}_t$ la force d'adhérence. L'introduction du frottement de Coulomb, avec $K_\mu$ le cône de frottement isotrope de Coulomb défini par :

$$ K_{\mu} = \{ (r_n, \mathbf{r}_t) \text{ tels que } f(r_n, \mathbf{r}_t) = \|\mathbf{r}_t\| - \mu r_n \le 0 \} \qquad (2) $$

conduit à un problème de complémentarité non linéaire. Dans la majorité des techniques numériques présentées dans la littérature, ces difficultés sont abordées au moyen d'approximations régularisantes. Les conditions de non interpénétrabilité des corps sont remplacées par des lois de répulsion suffisamment raides qui entrent en jeu lorsque deux corps s'approchent. De même, la loi de frottement sec de Coulomb est usuellement régularisée. Une façon plus rigoureuse de modéliser une telle loi de contact multivoque est l'utilisation de la théorie moderne des fonctions convexes et du calcul sous-différentiel. A partir de ce formalisme mathématique, J.-J. Moreau propose d'introduire le concept de pseudo-potentiel, qui est convexe et semi-défini positif [9]. La condition de Signorini et la loi de Coulomb s'expriment alors sous forme de loi de sous-normalité :

$$ -\dot{\mathbf{u}}_t \in \partial \Psi_{K_{r_n}}(\mathbf{r}_t) \iff \mathbf{r}_t \in \partial \varphi_{r_n}(-\dot{\mathbf{u}}_t). \qquad (3) $$

G de Saxcé [2] a proposé d'abandonner l'idée d'une séparation en deux pseudo-potentiels duaux et postulé l'existence d'une fonction unique des variables duales appelée « bipotentiel de contact » :

$$ b_c(-\dot{\mathbf{u}}, \mathbf{r}) = \Psi_{\mathbb{R}^-}(-\dot{u}_n) + \Psi_{K_\mu}(\mathbf{r}) + \mu t_n \| - \dot{\mathbf{u}}_t \| \qquad (4) $$

où $\Psi_{\mathbb{R}^-}$ désigne la fonction indicatrice de $\mathbb{R}^- = ]-\infty, 0]^1$ et $\Psi_{K_\mu}$ est une fonction indicatrice qui prend la valeur zéro si $\mathbf{r}_t \in K_\mu$ et $+\infty$ sinon. Elle doit être biconvexe et satisfaire une inégalité généralisant celle de Fenchel :

$$ \forall -\dot{\mathbf{u}}, \mathbf{r} \in \mathbb{R}^3, \Psi_{\mathbb{R}^-}(-\dot{u}_n) + \Psi_{K_\mu}(\mathbf{r}) + \mu r_n \| - \dot{\mathbf{u}}_t \| \ge -(\dot{\mathbf{u}}_t \mathbf{r}_t + \dot{u}_n r_n) \qquad (5) $$

D'autre part, il faut que les couples extrémaux du bipotentiel vérifient la loi de contact avec frottement sec, c'est-à-dire que le bipotentiel vérifie l'égalité :

$$ \mu r_n \| - \dot{\mathbf{u}}_t \| = -(\dot{\mathbf{u}}_t \mathbf{r}_t + \dot{u}_n r_n). \qquad (6) $$

<footnote n="1">¹Cette fonction prend la valeur zéro si $-\dot{\mathbf{u}}_n \in \mathbb{R}^-$ et $+\infty$ autrement.</footnote>
Contrairement à l'approche usuelle, la méthode du bipotentiel conduit à un seul principe variationnel et une seule inégalité. A partir de l'algorithme d'Uzawa, on aboutit à un algorithme de résolution de la loi constitutive fondé sur le schéma prédicteur-correcteur suivant :

$$ 
\begin{aligned}
\text{prédicteur :} \quad & \tau^{i+1} = \mathbf{r}^i - \rho^i [\dot{\mathbf{u}}_t^i + (\dot{\mathbf{u}}_n^i + \mu^i) \parallel \dot{\mathbf{u}}_t^i \parallel] \cdot \mathbf{n}, \\
\text{correcteur :} \quad & \mathbf{r}^{i+1} = \text{proj}(\tau^{i+1}, \mathbf{K}_\mu)
\end{aligned} 
\qquad (7) $$

où la projection sur le cône de Coulomb conduit, suivant la valeur de $\tau$, à l'un des trois statuts : non contact, contact avec adhérence ou contact avec glissement (voir figure 1).

<figure><img src="image_2.png" /><figcaption>Figure 1 – Projection de la prédiction sur le cône de Coulomb</figcaption></figure>

## 3 Lois de comportement homogénéisées

Alors que le tenseur des contraintes est depuis longtemps bien défini pour un milieu continu, la définition d'un tenseur des contraintes moyen à partir des forces s'exerçant à l'échelle du grain est beaucoup moins claire pour un milieu granulaire. Il existe plusieurs définitions différentes qui ne conduisent pas en général à un tenseur des contraintes moyen symétrique. De plus, il est important de noter que tous ces travaux reposent sur une analyse statique du problème, et ne sont valables que pour un arrangement quasi-statique de milieux granulaires.

Nous avons proposé une définition générale du tenseur des contraintes moyen pour un milieu granulaire, valable en statique comme en dynamique [1, 5, 3]. Cette définition prend en compte, non seulement les forces de contact, mais également les forces de volume (gravité et inertie) s'exerçant à l'échelle des grains. Le tenseur des contraintes moyen ainsi défini possède les propriétés d'un tenseur des contraintes de Cauchy : il est symétrique et invariant par translation.

Afin de définir le tenseur des contraintes dans un milieu granulaire, on procède comme suit. Cette définition est valable pour des grains quelconques, qui peuvent être déformables ou non. On suppose juste pour simplifier que les contacts entre grains sont ponctuels.

i) On se donne un Volume Elementaire Représentatif (VER) de volume $V$, contenant $N_p$ grains « p » (voir figure 2), tel que la frontière du VER ne coupe aucun grain.

ii) Chaque grain $V_p$ est supposé en équilibre dynamique sous l'action des forces de volume $\vec{f}$ ($\vec{f} = \rho(\vec{g} - \vec{x})$) incluant les forces d'inertie, et des forces de contact $\vec{r}_\alpha$ localisées aux points de contact $\vec{x}_{c\alpha}$ entre les grains. Ainsi, la contribution de chaque grain $V_p$ au tenseur des contraintes moyen est donnée par :

$$ \Sigma_p = \frac{1}{V_p} \left[ \int_{V_p} \rho (\vec{g} \otimes \vec{x} - \ddot{\vec{x}}) dV + \sum_{\alpha} \vec{r}_{\alpha} \otimes \vec{x}_{c\alpha} \right] \qquad (8) $$

où $\vec{g}$ désigne la gravité, et le point la dérivée par rapport au temps.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Le VER est un système de $N_p$ particules</figcaption></figure>

iii) Compte tenu de la propriété d'aditivité des moments internes, le tenseur des contraintes moyen $\Sigma$ sur le VER est ensuite calculé comme suit :

$$ \Sigma = \frac{1}{V} \sum_{p \in V} V_p \Sigma_p = \frac{1}{\sum_{p \in V} V_p (1+e)} \sum_{p \in V} V_p \Sigma_p \quad (9) $$

où $e$ désigne l'indice des vides.

Il convient de souligner l'importance des propriétés du tenseur des contraintes moyen défini par (8) pour le calcul sur un grand nombre de grains. Tout d'abord, la propriété d'additivité des moments internes permet de ramener le calcul du tenseur des contraintes moyen du VER à la contribution de chacun de ses grains. Ensuite, l'invariance par translation permet de choisir la meilleure origine pour calculer le tenseur des contraintes moyen de chacun des grains. Et ceci afin de réduire les erreurs numériques qui seraient dues à un choix de l'origine trop éloignée des grains. D'autre part, le choix de l'origine au centre de chaque grain permet d'annuler la contribution de la gravité.

Nous proposons de montrer maintenant, sur l'exemple analytique suivant, l'importance des effets dynamiques et en particulier des forces d'inertie, dans la contribution d'un grain au tenseur des contraintes moyen. Afin de rendre le calcul analytique possible, considérons un ensemble de grains (voir figure 3) s'écoulant sur un plan incliné et isolons, entre deux chocs consécutifs, un grain roulant sans glisser sur le plan incliné.

Soit $V_b$ un grain homogène de forme cylindrique roulant sans glisser sur un plan incliné d'un angle $\varphi$ (voir Fig. 3). En dimension 2, le cylindre de rayon $a$, de masse volumique $\rho$ et de masse totale $m_b = \rho V_b$, est assimilé à un disque. D'autre part, on considère un repère orthométrisé $R_p = (0, \vec{e}_1, \vec{e}_2, \vec{e}_3)$ tel que $\vec{e}_1$ et $\vec{e}_2$ soient parallèles au plan incliné. On considère également le repère $R_b = (c, \vec{e}'_1, \vec{e}'_2, \vec{e}'_3)$ attaché à la particule, animé de la vitesse angulaire $\vec{w} = \theta \vec{e}_3$ par rapport à $R_p$. Ainsi, la vitesse d'un point courant de la particule, de position $\vec{x}$ dans $R_p$, est égale à $\vec{x} = \vec{c} + \vec{w} \wedge (\vec{x} - \vec{c})$ où $\vec{c}$ désigne la position du centre de gravité de la particule.

En utilisant alors (8) pour calculer la contribution de la bille au tenseur des contraintes moyen, on a :

$$ \Sigma_b = \frac{1}{V_b} \left( \int_{V_b} \left[ \rho(\vec{g} - \vec{\dot{x}}) \right] \otimes \vec{x} dV + \vec{r} \otimes \vec{x}_c \right) = \Sigma_r + \Sigma_s \quad (10) $$

où $\Sigma_r$ désigne la contrainte due à la réaction de contact :
<figure><img src="image_4.png" /><figcaption>Figure 3 – Particule rigide roulant sans glisser sur un plan incliné.</figcaption></figure>

$$ \Sigma_r = \frac{1}{V_b} \vec{r} \otimes (\vec{x}_c \alpha - \vec{c}), \qquad (11) $$

et $\Sigma_s$ désigne la contrainte due aux efforts de rotations

$$ \Sigma_s = -\frac{1}{V_b}(I_c \cdot id - J) \left( (j(\vec{w}))^2 - j(\dot{\vec{w}}) \right) = \frac{\ddot{\theta}}{V_b} \begin{pmatrix} 0 & -J_{22} \\ J_{11} & 0 \end{pmatrix} - \frac{\dot{\theta}^2}{V_b} \begin{pmatrix} J_{11} & 0 \\ 0 & J_{22} \end{pmatrix} = \Sigma_1 + \Sigma_2, \quad (12) $$

où $J$ est le tenseur d'intertie de la particule, $I_c$ son moment d'inertie par rapport à l'axe $\vec{e}_3$, $id$ est l'identité de $\mathbb{R}^2$, et $j(\vec{w}) \cdot \vec{u} = \vec{w} \wedge \vec{u}$.

Dans le cas d'un cylindre $b$ roulant sur un plan incliné, en utilisant les inconnues du problème, $(\vec{c}, \theta)$, ainsi que les composantes $r_1$ et $r_2$ de la réaction de contact données par la résolution des équations de la dynamique, on peut montrer que l'on obtient (voir [1, 5] pour les calculs détaillés) :

$$ \vec{r} = r_1 \vec{e}_1 + r_2 \vec{e}_2 = m_b g \left( -\frac{1}{3} \sin\phi \vec{e}_1 + \cos\phi \vec{e}_2 \right), \quad \ddot{\theta} = \frac{2r_1}{m_b a}, \quad \ddot{c}_1 = \frac{2}{3} g \sin\phi. \qquad (13) $$

Dans ce cas, le tenseur des contraintes moyens d'un grain cylindrique roulant sans glisser sur un plan incliné est donné par

$$ \Sigma_b = \rho g a \left\{ \begin{pmatrix} 0 & \frac{\sin\phi}{6} \\ \frac{\sin\phi}{6} & -\cos\phi \end{pmatrix} + \frac{\mathcal{L}}{a} \begin{pmatrix} \frac{\sin\phi}{3} & 0 \\ 0 & \frac{\sin\phi}{3} \end{pmatrix} \right\}, \qquad (14) $$

où $\mathcal{L} = c_1(t) - c_1(0) = \frac{gt^2}{3} \sin\phi$.

On constate sur cet exemple, que la contribution $\Sigma_r$ due aux seules forces de contact, et qui correspond au terme classique de Weber, n'est pas symétrique. C'est la contribution $\Sigma_s$, due aux efforts d'inertie, qui est essentielle pour assurer la symétrie du tenseur des contraintes et obtenir un tenseur des contraintes moyen au sens de Cauchy.

D'autre part, $\mathcal{L}$ pouvant être identifié au libre parcours moyen entre deux chocs, on peut voir que le tenseur des contraintes moyen $\Sigma_b$ se décompose en deux termes. Dans le premier terme, les composantes $\Sigma_{12}$ sont positives et correspondent aux effets de friction lors du contact. La
composante Σ₂₂ est négative et correspond à une contrainte de compression, due à la réaction de contact. Le second terme de l'expression (14) de Σₚ, qui correspond au terme Σ₂ de (12), s'identifie à une pression hydrostatique positive qui correspond à un état de traction dû aux effets centrifuges. Dans le cas d'un matériau granulaire dense où le libre parcours moyen est petit, on a L << 1 et ce second terme est négligeable. Au contraire, dans le cas d'un matériau granulaire lâche, il ne peut pas être négligé. Afin d'illustrer cette propriété, la figure 4 montre un écoulement granulaire polydisperse de 480 particules de rayon 1.5mm(±0.5) sur une pente inclinée à 26 degrés avec un coefficient de frottement entre particules de μ = 0.5. La figure représente, pour chaque particule p le rapport de la norme de Frobenius de la matrice Σ₂ divisée par la norme de la contrainte Σₚ. Au plus la couleur est claire, au plus le rapport est grand. Nous observons que cette quantité est bien représentative des parties fluides et solides de l'écoulement.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Ecoulement dense (couleur foncée)/lache (couleur claire)</figcaption></figure>

# 4 Application

## 4.1 Cas test

Dans un premier temps, afin de caler le modèle numérique pour des contacts multiples, nous avons procédé à des essais sur des échantillons constitués de matériau analogique homogène. Soit un ensemble de cylindres rigides, homogènes confinés dans une boîte constituée de deux parois verticales et de deux parois horizontales immobiles. Chaque cylindre est sujet à l'accélération de gravité g et à des efforts de contact dus aux parois et aux cylindres voisins, Figure 5.

<figure><img src="image_6.png" /><figcaption>Figure 5 – Efforts normaux pour un arrangement (a) rectangulaire - (b) triangulaire - (c) boulets de canon</figcaption></figure>

L'algorithme modélise bien le phénomène de “descente des charges” caractéristique de l'arrangement rectangulaire. Pour un assemblage triangulaire, il met en valeur une caractéristique très importante dans les milieux discrets : la résultante des efforts de contact sur la paroi horizontale
inférieure n'est pas équivalente au poids total des cylindres mais inférieure, car les parois latérales subissent des efforts de contact normaux et tangentiels. Pour un assemblage dit "boulets de canon" on constate que c'est l'énergie potentielle de gravitation dont l'énergie potentielle est un minimum local qui assure la stabilité.

## 4.2 Stick-slip

Afin de mettre en évidence l'influence du frottement, nous avons choisis d'effectuer l'expérience numérique suivante : on considère un système de particules hétérogènes dans un récipient. Dans un premier temps, l'échantillon a été préparé en déposant les particules par gravité. Dans un second temps, nous imposons à la paroi du fond une vitesse de $5 \cdot 10^{-5} \text{ms}^{-1}$ vers le haut. Nous avons effectué une série de tests en modifiant le contact frottant particule-particule et particule-paroi. On constate pour des contacts non frottants que la force reste constante tout au long de la simulation. Par contre pour des contacts particules-parois frottant ($\mu = 0.8$), mais sans frottement entre particules, on constate que l'évolution de la force est très chaotique au début de la simulation et devient relativement régulière par la suite (figure (6a)). Si maintenant on a des contacts particule-particule et particule-parois frottants ($\mu = 0.8$), on constate qu'il y a une montée en charge de l'intensité de la force suivi d'une décharge (figure (6b)).

<figure><img src="image_7.png" /><figcaption>Figure 6 – Evolution de la force en fonction du temps : (a) contact particule-paroi frottants ($\mu = 0.8$), (b) contact particule-particule et particule-paroi frottants ($\mu = 0.8$)</figcaption></figure>

## 4.3 Polarisation des efforts de contacts

On se propose ici de déterminer si les efforts de contacts au niveau de la paroi sont dans les conditions dite " de Coulomb" à savoir sur la surface de glissement. Nous avons considéré un échantillons de 350 particules de rayon $1.5 \text{mm}(\pm 10)$. Le pas de temps utilisé est de $10^{-3}$. Le frottement statique est de 0.25. La première étape a consisté à sédimenter l'échantillon sous la seule force de gravité et à "relâcher" le système afin de se mettre dans les conditions de Coulomb, Figure 7a. A la fin du "relâchement", on calcule les contraintes moyennes au niveau des parois. On compare ensuite les résultats obtenus à ceux donnés par l'approche continue dite de Janssen, Figure 7b.

## 4.4 Perspectives

Cette étude préliminaire nous permettra de considérer par la suite un écoulement cisaillé où la quantification du phénomène de polarisation nous permettra à terme de proposer un modèle de couche limite pour ce type d'écoulement, à l'image de ce qui a été réalisé pour la modélisation des joints dans le cas de maçonneries [8]. Ce type d'étude a déjà été réalisé dans le cas de milieux
<figure><img src="image_8.png" /><figcaption>Figure 7 – (a) Polarisation des efforts de contact - (b) Comparaison Janssen - Eléments Discrets</figcaption></figure>

non cisaillé, comme par exemple dans un tambour tournant avec la méthode des éléments discrets
[10] ou pour des modélisations continues des tas de sable sec [4].

**Remerciements.** Cette étude a été mené au sein du projet » Grain De Sable », qui a été financé par l'Agence Nationale pour la Recherche (programme JC05-41831).

Références

[1] B. Chetouane. Approche combinée d'éléments finis/éléments discrets pour la modélisation des structures maçonnées, *thèse de doctorat*, Université de Montpellier 2, 2004.

[2] G. de Saxcé, Z. Q. Feng. New inequality and functional for contact with friction : The implicit standard material approach, *Mechanics of Structures and Machines*, **19**, 3, 301-325, 1991.

[3] G. de Saxcé, J. Fortin, O. Millet. About the numerical simulation of the dynamics of granular media and the definition of the mean stress tensor *Mechanics of Materials*, **36**, 1175–1184, 2004.

[4] S. Dumont, N. Igbida. On a Dual Formulation for the Growing Sandpile Problem, *European Journal of Applied Mathematics*, **20**, 169–185, 2009.

[5] J. Fortin. Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires, *thèse de doctorat*, Université de Lille 1, 2000.

[6] J. Fortin, O. Millet, G. de Saxcé. Numerical simulation of granular materials by an improved discrete element method, *Int. J. Numer. Meth. Engng*, **62**, 639–663, 2005.

[7] M. Jean. The Non Smooth Contact Dynamics Method, *Compt. M. Appl. Math. Engrg*, **177**, 235–257, 1999.

[8] F. Lebon, R. Rizzoni. Asymptotic study of soft thin layer : the non convex case. *Mechanics of Advanced Materials and Structures*, **15**, 12–20, 2008.

[9] J. J. Moreau. Sur les lois de frottement, de plasticité et de viscosité, *C. R. Acad. Sci.*, t. 271, Série 1, 608–611, 1970.

[10] M. Renouf, D. Bonamy, F. Dubois, P. Alart Numerical simulation of two-dimensional steady granular flows in rotating drums : on surface flow rheology, *Physics of fluids*, **17**, 13303–13315, 2005.