Le bruit de frottement entre deux surfaces rugueuses est dû à la vibration verticale engendrée par les impacts inter-aspérités de deux solides frottés. Les expérimentations ne donnent pas l?accès aux grandeurs mécaniques locales comme l?intensité de choc, la déformation d?aspérité, la pression locale.?Afin de révéler la relation entre les derniers et les paramètres macroscopiques, on propose une approche numérique 1D basée sur la méthode de décomposition modale, les formulations de Lagrange et pénalité, et certains schémas d'intégration

DANG, Viet Hung

LE BOT, Alain

I-Revues

20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
Simulation numerique directe en dynamique du contact
rugueux
Viet Hung DANGa, Alain LE BOT a
a. Laboratoire de Tribologie et Dynamique des Syste`mes, Ecole Centrale de Lyon, 36 Avenue Guy de
Collongue, 69134 Ecully Cedex
Re´sume´ :
Le bruit de frottement entre deux surfaces rugueuses est duˆ a` la vibration verticale engendre´e par
les impacts inter-aspe´rite´s de deux solides frotte´s. Les expe´rimentations ne donnent pas acce`s aux
grandeurs me´caniques locales comme l’intensite´ de choc, la de´formation d’aspe´rite´, la pression locale...
Afin de re´ve´ler la relation entre ces derniers et les parame`tres macroscopiques, on propose une approche
nume´rique 1D base´e sur la me´thode de de´composition modale, les formulations de Lagrange et pe´nalite´,
et plusieurs sche´mas d’inte´gration temporelle.
Abstract :
The origin of the roughness noise is the shock-induced normal vibration of antagonist asperities.
Several experiments demonstrate the logarithm law between the speed and the radiated acoustic power
with the arithmetic roughness of the surface. In order to analyze the relations between the dynamical
characteristics and the macro parametes, we propose a numerical approach with model 1D based on
the general equations of beam vibration. Two methods to determine the contact force are used either
the penalty or Lagrange multiplier with several integration schema
Mots clefs : simulation nume´rique ; contact dynamique ; Bruit de frottement
1 Introduction
Lorsque deux objets sont frotte´s l’un contre l’autre, ils engendrent du bruit de frottement. Ces bruits
se rencontrent tre`s souvent dans la vie quotidienne comme : le brossage, le bruit entre les mains, le
freinage de voiture, le crissement de porte. Du point de vue physique, on peut les classer en au moins
deux grandes familles selon les pressions de contact : les bruits d’instabilite´s me´caniques comme le
stick-slip avec de grandes pressions et le bruit de rugosite´ avec de faibles pressions de contact. Dans
ce dernier cas, le me´canisme ge´ne´rateur du bruit est intimement lie´ a` la rugosite´ des surfaces. Lors du
frottement des solides, les surfaces rugueuses engendrent des impacts inter-aspe´rite´s qui provoquent
des vibrations. Ces vibrations rayonnent ensuite le bruit dans l’environnement acoustique. Le proble`me
canonique que nous cherchons a` re´soudre nume´riquement est donc constitue´ de deux profils rugueux en
interaction plan-plan comme sur figure 1. Le profil du haut se de´place horizontalement avec une vitesse
Figure 1 – Glissement de deux profils rugueux. Les altitudes hi sont compte´es de manie`re oppose´e.
L’e´cart entre les deux profils est δ. A` l’instant t, le de´calage horizontal est impose´ ∆ = V t.
constante V, alors que celui du bas est fixe´ aux deux extre´mite´s. La forme des profils est donne´e par
1
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
l’altitude hi(x) de l’ensemble des points de la surface. L’e´cart initial entre les deux profils est δ. Dans
cet article, on pre´sente une approche nume´rique par l’utilisation de la technique de de´composition
modale pour re´soudre l’e´quation du mouvement, les me´thode de pe´nalisation et des multiplicateurs de
Lagrange pour la gestion du contact, et diffe´rents sche´mas d’inte´gration temporelle. Les calculs sont
effectue´s avec les hypothe`ses suivantes :
– Mode`le de poutre d’Euler-Bernoulli en flexion transverse seulement sans de´formation longitudinale.
– A chaque pas de temps, la position horizontale des points de la poutre supe´rieure est impose´e.
– Condition de non-pe´ne´tration (Kuhn-Tucker) des aspe´rite´s en contact.
2 Equation des poutres
L’e´quation du mouvement des poutres en flexion est,
Di∆
2ui +mi
∂2ui
∂t2
= f(t, x), (1)
ou` ui est la de´forme´e transverse de la poutre et f le champ de force de contact s’appliquant a` la
poutre. Di = EiIi est la rigidite´ de la poutre, Ei module d’Young et Ii moment d’inertie, mi est la
masse line´aire. Les conditions aux limites pour une poutre en appuis simples aux extre´mite´s x = 0 et
x = Li,
ui(0) = u
′′
i (0) = ui(Li) = u
′′
i (Li) = 0. (2)
Les modes ψki avec k = 0, 1, ... de la poutre i sont de´finis de sorte que,∫ Li
0
ψki (x)ψ
l
i(x)dx = δ
kl (3)
Di∆
2ψki = mi(ω
k
i )
2ψki (4)
avec
ψki (x) =
√
2
Li
sin(k + 1)
pix
Li
et ωki =
√
Di
mi
(
(k + 1)pi
Li
)2
(5)
A` chaque instant t, on peut de´velopper la de´forme´e ui dans la base modale,
ui(t, x) =
∞∑
k=0
Uki (t)ψ
k
i (x) (6)
Il s’agit de trouver et re´soudre les e´quations sur les amplitudes modales Uki .
mi
[
U¨ki + (ω
k
i )
2Uki
]
= F ki (t) (7)
avec F ki la force modale de´finie par F
k
i (t) =
∫ Li
0
f(t, x)ψki (x)dx.
Dans le cas d’une structure amortie, l’e´quation du mouvement dans la base modale devient
mi
[
U¨ki + 2ζ
k
i ω
k
i U˙
k
i + (ω
k
i )
2Uki
]
= F ki (t) (8)
Nous e´crivons l’e´quation du mouvement sous la forme d’un syste`me d’e´quations diffe´rentielles du pre-
mier ordre.
[
U˙ki
V˙ ki
]
=
[
0 1
−
(
ωki
)
2 −2ωki ζ
k
i
] [
Uki
V ki
]
+
[
0
F ki
]
(9)
En posant y =
[
U(t)
V (t)
]
ces e´quations se retrouvent la forme ge´ne´rale du proble`me de Cauchy :
y˙ = g(y) ; y(0) = 0 (10)
Afin de re´soudre cette e´quation diffe´rentielle, on discre´tise en temps et en espace. Le pas de temps est
τ et le vecteur temps est tm = mτ . Le pas de distance est χi et le vecteur position est x
n
i = nχi. Les
de´forme´es, force et distance deviennent um,ni = ui(t
m, xni ), f
m,n = f(tm, xni ) et e
m,n = e(tm, xni ). Les
quantite´s modales deviennent Um,ki = U
k
i (t
m) et ψk,ni = ψ
k
i (x
n
i ). Enfin, les profils sont h
n
i = hi(x
n
i ).
2
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
3 Calcul de la force de contact
On conside`re un point quelconque sur le profil supe´rieur, sa distance normale avec le deuxie`me profile
est donne´ par la formule
g(t, x) = h1(x) + u1(t, x) + h2(x−∆) + u2(t, x−∆)− δ (11)
Si f de´signe la pression de contact en ce point, la condition de non-pe´ne´tration se manifeste donc sous
forme suivant :
g ≥ 0 ; f ≤ 0 ; g.f = 0 (12)
S’il n’y a pas de contact, g > 0, la pression de contact est naturellement nulle. Au contraire, si le
contact se produit, il engendre une pression qui a tendance a` e´carter les deux corps solides.
Pour la simulation nume´rique du proble`me de contact, les deux me´thodes les plus populaires sont celle
de pe´nalite´ et celle des multiplicateurs de Lagrange.
La me´thode de pe´nalite´ : dans cette me´thode, on autorise une le´ge`re pe´ne´tration des points de contact
dans le corps solide antagoniste, qui vaut e=H(g).ou` H(g) = 0 si g < 0 et H(g) = g si g > 0. la pression
de contact est fournie par une relation line´aire : f(t, x) = κe(t, x) ou` κ est le coefficient de pe´nalite´.
Sa valeur obtenue empiriquement ne soit pas eˆtre trop faible en raison de pe´ne´tration irre´aliste et ne
doit pas eˆtre trop grande ce qui va conduire une instabilite´ nume´rique.
La me´thode des multiplicateurs : on respecte parfaitement le condition de non-pe´ne´tration, g ≥ 0,
juste avant l’instant de contact, on introduit λ qui s’appelle un multiplicateur de Lagrange. Ce mul-
tiplicateur jouant le roˆle de pression de contact fictive va empeˆcher les de´formations conduisant a` la
pe´ne´tration.
L’avantage de la me´thode pe´nalite´ est sa facilite´ de mise en oeuvre, mais il faut choisir de´licatement
la valeur du coefficient et accepter une pe´ne´tration non-physique. Alors que la me´thode des multi-
plicateurs de Lagrange, donne des solutions tre`s proches de re´alite´, mais elle introduit des variables
supple´mentaires qui alourdissent conside´rablement le calcul.
4 Inte´gration temporelle
Me´thode Diffe´rences Euler Euler RK2 RK4 RK8
centre´es modifie´
Pas τ < 2
ωk
i
instable τ < 2
ωk
i
τ <
√
2
ωk
i
du temps
Erreur o (τ) o (τ) o
(
τ2
)
o
(
τ4
)
o
(
τ5
)
locale
Proprie´te´ Explicite Explicite Explicite Explicite Explicite Explicite
pas mul-
tiple
un pas un pas un pas un pas un pas
Table 1 – Sche´ma d’inte´gration temporelle
Comme le proble`me se manifeste sous la forme d’un proble`me ge´ne´ral de Cauchy, on va utiliser
les sche´mas inte´grations temporelles ordinaires pour le re´soudre nume´riquement. Ces sche´mas nous
permettent d’obtenir la de´formation, la vitesse, la force des points sur le profil a` chaque pas de temps.
Dans cette e´tude, on propose 6 sche´mas d’inte´gration temporelle. L’e´valuation d’un sche´ma se base
sur 3 crite`res : la consistance, la stabilite´ et la convergence. Le sche´ma est consistant si la somme de
l’erreur locale tend vers 0 lorsque le pas de temps devient petit. D’autre part, si la propagation des
erreurs d’arrondi peut eˆtre controˆle´e au fil des calculs, le sche´ma est stable. Un sche´ma consistant et
stable est dit convergent c’est-a` -dire que l’erreur globale tend vers 0 lorsque le pas devient petit. Le
tableau 1 re´sume le nom et les caracte`res des sche´mas utilise´s. En se basant sur la clarte´ des proprie´te´s
et de la facilite´ d’implantation, on a choisi le sche´ma Euler modifie´ comme le sche´ma principal pour
effectuer les calculs.
3
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
5 Validation du programme Ra1D
Un cas test tre`s simple entre deux surfaces glisse´es a` vitesse constante (voir figure 2) est utilise´ pour
comparer les deux programme ra1D (code maison) et Abaqus explicit. La solide infe´rieur est appuye´
aux deux extre´mite´s. Le maillage utilise´ est de type triangulaire quadratique, 6 noeuds par e´le´ment
pour abaqus. Il a e´te´ raffine´ au niveau des aspe´rite´s des surfaces en contact. Les hauteurs d’aspe´rite´s
sont de l’ordre de 1e−6 (m). Pas spatial au niveau de l’interface χ = 1e−4m.
Les dimensions de deuxs solides sont L = 0.01 m (longueur), H = 0.0005 m (hauteur). Les ca-
(a) Maillage du mode`le (b) Maillage agrandi
Figure 2 – Maillage du cas test sixChocs
racte´ristiques me´caniques du mate´riau sont E = 69000Mpa pour le module d’Young, µ = 0.33 pour le
coefficient de Poisson, et ρ = 2000kg/m3 pour la masse volumique. Coefficient de pe´nalite´ : κ = 1e13Pa.
Pas du temps τ = 2e−9s et dure´e de simulation T = 0.001s. Dans ce cas test, on ne´glige la force de
gravite´, la friction, et l’amortissement structural. La figure 3 montre le de´placement et la pression de
contact au point M le plus haut du profil supe´rieur. Une bonne cohe´rence est observe´e entre Ra1D et
Abaqus. Les temps CPU, il sont 280s avec Abaqus et seulement 30s avec Ra1D.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−1
−0.5
0
0.5
1
x 10−7
Temps (s)
D
ep
la
ce
m
en
t (
m
)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−1
−0.5
0
0.5
1
x 10−7
Temps (s)
D
ep
la
ce
m
en
t (
m
)
(a) De´placement d’un point
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−6
−4
−2
0
x 105
Temps (s)
Pr
es
sio
n 
de
 co
nt
ac
t (N
)
 
 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x 10−3
−10
−8
−6
−4
−2
0
x 105
Temps (s)
Pr
es
sio
n 
de
 co
nt
ac
t (N
)
 
 
Ra1D
Abaqus
(b) Pression de contact en un point
Figure 3 – Comparaison de re´sultat entre Ra1D avec Abaqus
Pour le deuxie`me cas test, on compare Ra1D avec une solution analytique dans le cas de la vibration
d’une poutre provoque´ par le de´placement d’une masse. La solution analytique est pre´sente´e dans
M.Olsson [5] Sur la figure 4, on trace la de´flection du milieu de la poutre au cours du temps avec
diffe´rentes vitesses de de´placement de la masse. Avec P est le poids de la masse mobile se de´plac¸ant a`
vitesse V, α = v T1
L
et T1 : pe´riode du premier mode de la poutre. Une concordance parfaite avec les
re´sultats obtenues par la solution analytique, a e´te´ observe´e.
6 Re´sultats de calculs
Les profils sont gaussiens et ge´ne´re´s nume´riquement. Les parame`tres de topographie sont indique´s
tableau 2. Les poutres sont en acier avec un module d’Young E = 210 GPa et une masse volumique
ρ = 7800 kg/m3. Elles ont pour longueur L = 120mm, e´paisseur b = 4mm. Le taux d’amortissement
modal est ξ = 0, 01. Les valeurs des fre´quences propres sont pre´sente´es tableau 3.
Les donne´es ne´cessaires pour la simulation sont pas de temps τ = 1e−8s, pas spatial χ = 8µm, dure´e
de simulation T = 1s.
Le premier re´sultat obtenu est que le temps CPU de Ra1D est beaucoup plus petit que celui d’Abaqus
(10h contre 7 jours). Sur la figure 5, un exemple d’audiogramme simule´ et sa densite´ spectrale de
puissance sont montre´s. On trouve que les pics les plus hauts correspondent aux fre´quences propres
4
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
−4.5
−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5
0
0.5
x 10−3
Position relative de la masse sur la poutre
D
ep
la
ce
m
en
t (m
)
 
 
a=0.125, ra1D
a=0.25,   ra1D
a=0.5,     ra1D
a=1.0,     ra1D
a=0.125, analytique
a=0.25,   analytique
a=0.5,     analytique
a=1.0,     analytique
(a) De´placement du point au milieu de la poutre
selon le temps.
(b) Masse de´place´e sur une poutre
Figure 4 – Comparaison Ra1D / solution analytique pour une masse mobile sur poutre
Rq
(µm)
Rt
(µm)
Rz
(µm)
Ra
(µm)
Rsk Rku RSm
(µm)
02 12.8 6.0 1.7 0.012 3.6 1574
04 27.8 12.2 3.5 -0.009 4.6 1520
08 53.0 23.5 6.8 0.041 3.6 1560
12 86.7 36.1 10.4 -0.019 3.8 1552
20 138 61 17.6 -0.125 4.2 1546
Table 2 – Proprietes topographiques
Mode 1 2 3 4 5 6
f(Hz) 1.31e+03 5.23e+03 1.18e+04 2.09e+04 3.27e+04 4.71e+04
Table 3 – Fre´quences propres des profils
Figure 5 – Audiogramme simule´
des poutres de´couple´es. Le niveau de la vibration est de´fini comme
Lv(dB) = 20log(
VRMS
V0
) (13)
VRMS =
1
Ntemps
∗
1
Nnoeud
ΣVi,j (14)
Ou` V0 = 1e
−9m/s est la vitesse vibratoire de re´fe´rence. VRMS est la vitesse moyenne quadratique, Vi,j
est la vitesse du noeud i a` l’instant j. N et Nt sont respectivement le nombre de noeuds et de pas de
temps. La relation entre le niveau de vibration et la rugosite´ ou la vitesse de glissement est montre´
dans fig 6 on trouve que Lv est proportionnel au log de la rugosite´ de la surface et e´galement de la
5
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
−1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4
110
120
130
140
150
160
170
Rapport entre Lv et logarithme de vitesse
Lv
(dB
)
Log(v)
 
 
Rq02, n= 1.12
Rq04, n= 1.11
Rq08, n=1.13
Rq12, n=1.05
Rq20, n=1.11
(a) Lv=f(V)
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
115
120
125
130
135
140
145
150
155
160
165
Lv
 (d
B)
log(Rq) (um) 
 
 
v=3cm/s
v=15cm/s
v=20cm/s
v=25cm/s
v=40cm/s
(b) Lv=f(Rq)
Figure 6 – Variation du niveau de la pression acoustique Lv(dB) en fonction de la vitesse et de la
rugosite
vitesse de frottement qui s’accorde avec la formule suivante :
∆Lv(dB) = 20 log
[(
V2
V1
)m]
avec m ≃ 1, 1 (15)
Ces re´sultats montrent une bonne cohe´rence avec les expe´riences re´alise´es dans H. Ben Abdelounis[4]
7 Conclusion
Le syste`me mode`le e´tudie´ ici permet de calculer pre´cise´ment la dynamique locale de deux surfaces
rugueuses en glissement relatif. L’e´quation du mouvement est re´solue nume´riquement graˆce a` la tech-
nique de de´composition modale et la double discre´tisation en temps et en espace. Plusieurs sche´mas
d’inte´gration temporelle ont e´te´ programme´s et compare´s. Le code obtenu, limite´ au cas plan, est plus
rapide que la me´thode des e´le´ments finis standard. Le logiciel est valide´ en comparaison avec Abaqus,
la formule analytique et les re´sultats expe´rimentaux
Re´fe´rences
[1] Le Bot A., Bou-Chakra E. 2009 Measurement of friction noise versus contact area of rough surfaces
weakly loaded. Tribology Letters. vol. 37, 273-281
[2] Le Bot A., Bou-Chakra E., Michon G. 2011 Dissipation of Vibration in Rough Contact. Tribology
letters vol. 41, 47-53
[3] BAILLET Laurent, SASSI Taoufik 2003 Numerical implementation of different finite element
methods for contact problems with friction. Comptes rendus. Me´canique vol. 331, 789-796
[4] Ben Abdelounis H., Le Bot A., Perret-Liaudet J., Zahouani H. 2010 The roughness effect on the
frequency of frictional sound. Wear A. vol. 268, 335-345
[5] OLsson M. 1991 On the fundamental moving load problem. Jounal of Sound and Vibration vol.
145, 299-307
[6] AKAY Adnan 2002 Acoustics of friction . The Journal of the Acoustical Society of America vol.
111, 1525-1548
[7] FORD I. J. 1993 roughness effect on friction for multi-asperity contact between surface. Journal
of physics. D, Applied physics vol. 26, 2219-2225
[8] STOIMENOV Boyko L., MARUYAMA Suguru, ADACHI Koshi, KATO Koji 2007 The roughness
effect on the frequency of frictional sound. Tribology international vol. 40, 659-664
6


Simulation numérique directe en dynamique du contact rugueux

http://documents.irevues.inist.fr/bitstream/2042/46316/1/cfm2011_538.pdf

Simulation numérique directe en dynamique du contact rugueux

Abstract

Similar works

Full text

Available Versions

I-Revues