Modélisation par éléments distincts

d’ouvrages en génie civil.

La méthode « Non Smooth Contact
Dynamics »

Robert Péralès* —Frédéric Dubois**—Marc Vinches*—Claude

Bohatier**

*
Ecole des Mines d’Ales. CMGD
6 Avenue de Clavières, 30319 Ales Cedex.
robert.perales@ema.fr, marc.vinches@ema.fr

**
Laboratoire de Mécanique et Génie Civil. Université de Montpellier 2. cc048
Place Eugène Bataillon, 34095 Montpellier cedex 05.
dubois@lmgc.univ-montp2.fr, bohatier@lmgc.univ-montp2.fr

RÉSUMÉ . Après une présentation des bases théoriques et bibliographiques relatives à la

méthode Non Smooth Contact Dynamics, trois structures réelles de génie civil sont
présentées. Elles ont la particularité d’illustrer le comportement d’ouvrages en pierre post-
contrainte, et d’un ouvrage antique, comportant plusieurs milliers de blocs, en en respectant
l’appareil.
ABSTRACT. After a recollection of the theoretical basis and main references concerning the
Non Smooth Contact Dynamics method, three real civil engineering structures are presented.
They have the peculiar features of a post-stressing use, and the modelling of a Roman bridge,
with thousands of tri-dimensional blocks, representing the original stone course.
MOTS-CLÉS : maçonnerie en grand appareil, précontrainte, mécanique non régulière, génie
civil.
KEYWORDS: large stone course masonry, pre-stressing, non smooth contact dynamics, civil
engineering.
 2

1. introduction

Avec le développement des programmes de préservation du patrimoine, un

intérêt accru est porté à l’étude des ouvrages maçonnés et à leur modélisation sous
différentes sollicitations. Dans ce domaine la structure peut être considérée comme
une collection de corps rigides ou déformables entre lesquels des relations
d’interaction, usuellement du contact avec frottement, sont susceptibles de s’établir
ou de se rompre. Ces types de milieux « divisés » composés d’un grand nombre de
corps peuvent être modélisés par les méthodes par éléments distincts (MED).
Le présent travail est le résultat d’une collaboration entre l’Ecole des Mines
d’Alès et le Laboratoire de Mécanique et Génie Civil de l’Université de Montpellier
II. Il s’appuie principalement sur le savoir faire développé par Michel Jean et Jean
Jacques Moreau dans le domaine de la simulation des grandes collections de corps
rigides ou déformables en contact.
De nombreux ouvrages de génie civil sont des structures maçonnées. Dans de
nombreux cas, ces structures sont étudiées comme des milieux continus. Cette
hypothèse n’est pas réaliste dans certaines applications comme la construction en
pierre massive à joints vifs, ou la construction en appareillages de pierres avec
mortier, de propriétés mécaniques faibles. Dans le premier cas de figure, on peut
classer des monuments antiques tels que l’aqueduc d’Arles, au franchissement du
vallon des Arcs à Fontvieille, ou des parties du pont du Gard, ou encore des arènes
de Nîmes.
Les figures suivantes présentent divers éléments de structure de ces monuments,
illustrant le fait qu’il apparaît réaliste de les modéliser partiellement, comme des
collections de corps discrets ou blocs.

Figure 1. Vue générale de l’aqueduc d’Arles et de l’aqueduc des moulins, au vallon des
Arcs, à Fontvieille.

Les modes de construction (faces de blocs planes, blocs montés à sec) ainsi que les
faciès de ruine de cet édifice ne permettent pas de considérer l’ouvrage comme
continu.
A la suite de la construction d’un deuxième aqueduc, à l’Est, dit aqueduc de
Barbegal, alimentant une meunerie industrielle, il semble que le premier aqueduc
d’Arles, ait été ruiné, et qu’un deuxième aqueduc ait été reconstruit, sur les piles du
 3

premier, en utilisant une technique de petit appareillage en façade, et utilisation

d’une maçonnerie de blocage à l’intérieur de la structure supérieure de l’ouvrage,
schématisé figure .

(10)

(9) (1)

(2)

(3)

(8)
(4)

(5) (6)

(7)

Figure 2. Modèle des maçonneries en grand appareil en partie inférieure, et en petit appareil
et maçonnerie de blocage en partie supérieure, du deuxième aqueduc d’Arles, au vallon des
Arcs, à Fontvieille.
La modélisation de tels ensembles par des modèles continus ne semble pas réaliste,
notamment pour ce qui concerne les parties de structure en grand appareil à joints
secs.
Dans le cadre des études entreprises sous l’impulsion de la Direction Régionale
des Affaires Culturelles de la région Languedoc-Roussillon, concernant les arènes de
Nîmes (fig. 3), il est apparu nécessaire de modéliser le bâtiment, de façon à obtenir
des plans de coupe, de recollement d’interventions et des indications réalistes en
termes de masses de maçonneries mises en œuvre, descentes de charges, suivis des
pathologies, etc. L’Ecole des Mines d’Alès a entrepris de réaliser une modélisation
géométrique réaliste du bâtiment, (fig. 4) à partir de plans, de relevés sur le terrain,
et d’une synthèse d’architecture réalisée par Myriam Fincker. Cette étape a donné
lieu à la simulation de descentes de charges, sur certaines géométries particulières,
en utilisant un modèle éléments discrets 3D.

Figure 3. Vue partielle des arènes de Figure 4. Arènes de Nîmes. Exemple de

Nîmes, détail de maçonneries en grand géométrie modélisée (CMGD, Ecole des
appareil. Mines d’Alès)
La relative complexité de la géométrie, réellement tridimensionnelle de l’édifice,
l’agencement particulier des blocs de pierre utilisés pour la partie la plus externe de
 4

l’édifice (galeries extérieures, et de l’attique), les techniques de construction des arcs

de décharge en façade, et des linteaux de la galerie extérieure du premier étage ne
sont évidemment pas prises en compte dans un modèle continu élastique. Nous
exposerons comment des arrangements de polyèdres tridimensionnels peuvent être
modélisés de façon plus réaliste par les méthodes aux éléments distincts (MED).
Nous allons rappeler les grandes lignes des méthodes par éléments distincts, et plus
particulièrement de la méthode NSCD. Nous présenterons quelques résultats
numériques sur trois structure réelles. Ces résultats on été obtenus avec le code de
calcul LMGC90.

2. Méthodes par éléments distincts: la méthode NSCD

Les méthodes par éléments distincts (MED) ont pour objet de modéliser un
milieu divisé au travers du comportement collectif de ses composants, chaque
composant pouvant par ailleurs être un objet complexe. Ces méthodes sont très
utilisées, par exemple, pour la modélisation des milieux granulaires (sable, roches,
grains dans des silos, etc.), des milieux condensés (gels), des milieux vivants
(cellules osseuses), certains milieux continus endommageables ou micro fracturés,
mais aussi les mécanismes (robotique).
2.1. Description d’un mouvement d’un corps en présence de « contact »

Les équations de la dynamique régissant le mouvement d’un corps en présence de

« contact » (mais en l’absence de choc) peuvent s’écrire sous la forme synthétique
suivante :

 M q&& = F (q, q& ) + P(t ) + r


conditions initiales et aux limites

où q désigne les variables de Lagrange (par exemple translations et rotations du

centre de gravité d’un corps rigide, ou déplacements d’un nœud d’un maillage), r
désigne la contribution des efforts de « contact » à cette équation (torseur ou réaction
nodale) et F ( q, q& ) + P (t ) représente la somme des efforts internes et externes.
Pour un corps déformable l’expression de l’équation de la dynamique peut être
obtenue par discrétisation par éléments finis du problème continu. Pour un corps
rigide cette équation signifie que le torseur des efforts extérieurs vaut le torseur
dynamique. On peut lui préférer d’autres formes comme les équations d’Euler
Newton :

 M u&& = P(t ) + r

I ω& + ω ∧ ( I .ω) = Mp(t ) + Mr
 5

Où u désigne le déplacement, ω la vitesse de rotation autour du centre d’inertie du

corps, P(t) et Mp(t) la résultante et le moment résultant des forces appliquées, r et
Mr la résultante et le moment résultant des efforts de contact, et M,I les matrices de
masse et d’inertie. La seconde équation de ce système est éventuellement non
linéaire à cause du second terme de gauche. Ce terme disparaîtra dans une
formulation 2D ou lorsque le corps est « isotrope » géométriquement. Il existe
d’autres formes d’écriture des équations du mouvement des solides rigides, par
exemple les équations de Lagrange, les formulations variationnelles (voir Brogliato
et al. (2002)).
2.2. Intégrateurs

Contrairement à la discrétisation en espace, la discrétisation en temps est perturbée

par l’introduction de conditions de « contact ». En effet en présence de chocs, ou ne
serait ce qu’en présence de contacts simultanés avec un découpage temporel
arbitraire, les dérivées précédentes doivent se comprendre au sens des distributions
(Moreau & Jean (1992)). On parlera dans ce cas de problème non régulier. Quel
que soit le schéma d’intégration nous cherchons à obtenir une solution discrète des
équations différentielles, i.e. à des instants : t0 < … < tn..
Dans le cas non régulier deux approches sont envisageables :

L’approche événementielle (event-driven) où le découpage temporel est
ajusté sur les évènements non réguliers. On essaie autant que possible de
séparer la partie régulière du mouvement (vol libre, contact établi) de la
partie non régulière du mouvement (choc, décollement). Cette méthode
permet d’utiliser quasiment n’importe quel intégrateur d’ordre élevé
(Abadie (1998)) mais nécessite de reformuler les conditions de contact en
accélération ce qui est contestable.
• L’approche pas à pas (time-stepping) où le découpage en temps est fixé de
manière arbitraire et non nécessairement constant. La méthode est capable
d’intégrer les équations du mouvement au cours d’un pas de temps donné.
Dans ce cas, l’ordre de l’intégrateur est nécessairement faible. La méthode
de Time-stepping est la mieux adaptée pour aborder le problème de
discontinuités temporelles.

3. Détection des intéractions

Le principe général utilisé pour la détection est le suivant:

détection d’une proximité grossière des contacteurs (pré détection). Partant
d’informations géométriques grossières, de valeurs de distance d’alerte on va
identifier les zones potentielles de contact.

test sur la visibilité des contacteurs. Une table de visibilité permet de faire
un trie dans les contacts potentiels en décidant arbitrairement que certains ne
 6

sont pas valables. De plus cela permet d’affecter une loi de contact
dépendant de la « couleur » des contacteurs en interaction (propriété
mutuelle de deux contacteurs).

détection précise du lieu de contact: point, ligne, face (détection locale).

calcul des informations locales: repère, interstice, vitesse relative,etc.
Dans toutes ces opérations c’est la détection grossière et le calcul précis du point
de contact qui s’avère le plus coûteux. On trouve un certain nombre de méthodes de
détection (pré ou locale):

méthode des boîtes de taille uniforme (pré). En se callant sur la taille des
particules et la valeur des distances d’alerte, … on vient faire un maillage de
l’espace. Grâce à ça on ne détectera que les particules situées dans une boîte
et ses voisines,

boites

Figure 5. Création du partitionnement en boîtes sur un échantillon.

Méthodes de Shadow-overlap (JJ Moreau), plan séparateur (Cundall), GJK.
Il s’agit de méthodes plus « mathématiques » de calcul de points de contact.
Elle se base sur des techniques de projection. Elles ont l’avantage d’être très
puissantes et assez génériques.
3.1. Cinématique et comportement des interactions

• Les lois de comportement au « contact » entre un contacteur candidat (C) et

un corps antagoniste (A) sont décrites par une relation point à point. On suppose
qu’on peut définir à chaque instant les deux points les plus proches des contacteurs
(A et C), et un repère local (n,t,s) porté par le contacteur antagoniste; la normale
étant dirigée vers le contacteur candidat ; (t,s) définit le plan tangent au contacteur
antagoniste. Les lois d’interaction les plus utilisées sont :

Pour le contact : loi de Signorini.
Dans sa forme classique la loi de Signorini relie
l’interstice à la réaction normale par une
RN relation de complémentarité:

g ≥ 0 , R N ≥ 0 , gR N = 0
UN,g
Pour le contact entre corps rigides JJ Moreau a montré
qu’en cas de contact « établi » cette relation était
 7

équivalente à une relation de Signorini exprimée en vitesse:

U N ≥ 0 , RN ≥ 0 ,U N RN = 0
Dans certains cas on lui préfèrera une version régularisée.

Pour le frottement : loi de Coulomb

RT
•
µRN Cette expression relie la vitesse
relative tangente à la réaction tangente, il s’agit
UT d’une loi à seuil :

 RT < µRN ⇒ UT = 0
− µRN 
RT ≤ µRN  R
RT = µRN ⇒ UT = −λ T ; λ ≥ 0
Coulomb
 RT

De même on peut lui préférer une loi régularisée, ou une version avec un coefficient
de frottement variable.
Il existe beaucoup d’autres lois, on pourra par exemple consulter à se sujet Cambou
et Jean (2001).
4. Stratégies et méthodes de résolution

On pourra trouver dans Brogliato et al. (2002) un excellent panorama sur la

simulation numérique des systèmes dynamiques non réguliers, et dans Jean (1993)
une description de plusieurs stratégies et méthodes de résolution.
Nous allons voir ici les deux derniers ingrédients entrant dans la composition de
nos « solveurs » : la façon de poser le problème de « contact » et la méthode
employée pour le résoudre. D’un point de vue formel, sur la figure 6, nous
pouvons schématiser le problème à résoudre de la façon suivante :

{q&} ← équation de la dynamique → {r} niveau global

∗
H (q) ↓ ↑ H (q)
{U }
α
← loi de contact → {R }
α
niveau local

Figure 6. Relations entre niveau local et niveau global

Au niveau global nous avons l’ensemble des systèmes d’équations de la dynamique
vérifiés par les corps. Au niveau local nous avons les lois de contact qui mettent en
relation des paires de nœuds portés par des contacteurs. Pour passer d’un niveau à
l’autre nous avons les matrices H et H* qui dépendent de la forme des zones de
contact, de la discrétisation, etc.
 8

H ∗ (q)
n n
(U , R) (U , R)
(q& , r ) t
H (q)
(q&i , ri )

a - cas rigide b- cas déformable

Figure 7. relation entres variables globales et locales.

Les approches de résolution seront donc les suivantes :

soit on cherche à résoudre le problème simultanément sur le niveau global et
local. C’est l’approche employée par Alart & Curnier (1991),

soit on résout au niveau global. C’est le principe des stratégies d’élimination
ou de substitution,

soit on résout au niveau local. C’est le principe des stratégies de
condensation.

Une fois l’approche de résolution choisie il reste à déterminer la méthode de
résolution.
4.1. Approche globale, Approche locale

Le principe de l’approche globale consiste à résoudre le système linéaire composé de

l’ensemble des inconnues globales tout en tenant compte des relations de contact.
• Partant de l’ensemble des systèmes d’équations de la dynamique des corps,
discrétisés en espace et en temps, on peut assembler le système global:

IM ( q& + − q& − ) = h r free + h r

• La matrice IM est diagonale par bloc, chaque bloc ayant la taille du système
d’équations définissant le mouvement de chaque corps. rfree représente le
vecteur résidu libre, obtenu sans tenir compte de la contribution du contact,
h est le pas de temps.
• A ce système d’équation vient se rajouter l’ensemble des relations de
contact. L’obtention d’un système ne dépendant plus des variables locales
est décrite dans Jean (1993) et ne sera donc pas reprise en détail ici. Le
principe général est globalement le suivant. En tenant compte des relations
cinématiques (matrices H) on élimine les variables cinématiques locales (g,
UT) pour ne conserver que les variables cinématiques globales (q, q& ).
Ensuite on utilise le système global pour éliminer r de l’expression
précédemment obtenue, ce qui nous donne un système d’inéquations ne
 9

dépendant que de q& . Pratiquement ces inéquations peuvent être résolues

avec n’importe quelle méthode de minimisation sous contrainte.
• D’un point de vue numérique, cette stratégie de résolution va permettre de
construire différentes méthodes avec résolution directe ou itérative.
• Dans l’approche locale, on va venir condenser (explicitement ou
numériquement) les équations de la dynamique au niveau des contacts. De
ce fait la résolution va se faire en deux temps tout d’abord au niveau des
relations de « contact » et ensuite une fois que toutes les relations sont
vérifiées au niveau des corps. Pour la suite on suppose qu’on est capable de
condenser le système global sous la forme suivante :

H∗ IM−1H h R − H∗ q& + = − H∗ ( IM−1hrfree + q& − ) ⇒ W hR − U = −U free

• On peut aussi l’écrire sous forme d’équation :

U α = U αfree ∑ wαβ h R β
nc
+
1

• Comme précédemment on va pouvoir distinguer des méthodes conduisant à

une résolution directe ou itérative. Pour avoir accès à certaines méthodes de
minimisation sous contrainte il est nécessaire de reformuler le problème
condensé (Jean (1993)) en y injectant les relations de « contact » ce qui va
nous donner accès à un système d’inéquations ne dépendant que des
réactions locales R.
5. Exemples de structures

On présente dans cette partie, trois structures réelles modélisées avec le code de
calcul LMGC90. Les blocs sont considérés comme rigides.
5.1. Escalier Ridolfi

Les caractéristiques de cet escalier, nous ont été données par la faculté d'architecture
de Bari. L'ouvrage, présenté à la 40° foire internationale de Vérone en 2005 sur les
techniques nouvelles appliquées à la pierre, est monté à sec et armé avec des câbles
en acier afin d’obtenir une précontrainte qui permette l’équilibre de la structure. Il
veut contribuer à montrer la résolution de certains problèmes qui concernent
l’emploi de la pierre dans l’architecture contemporaine.
• Contraintes et exigences
L’escalier est constituée de 20 marches. La loi de comportement, au niveau de
l’interface entre les marches, est de type frottement sec de Coulomb, standard, avec
un coefficient de frottement de 0.8.(A.Baudoin & al, 2006)
 10

Pour simuler une précontrainte, on a inséré un ressort entre les marches (celui-ci
modélise une raideur active en traction et inactive en compression).
Enfin, la première marche est encastrée sur un socle (la liaison est de type surface –
surface, avec une très forte cohésion).
Pour une raideur k du ressort de l’ordre de 6,12MN/m, une pré-déformation de 0.05
et une force de pré-contrainte équivalente de l’ordre de 107kN, on obtient après
5000 itérations (pas de temps de 10-4s) :

Figure 8. Représentation de la répartition Figure 9. Représentation de la répartition

des efforts des vitesses des marches
La répartition des efforts au niveau de chaque marche est homogène, tous les blocs
sont contraints de la même façon. Les vitesses des blocs sont de l'ordre de 10-5 m/s
donc négligeables : la structure est stabilisée.
La prise en compte d’une précontrainte trop faible, ne permet pas d’assurer un
équilibre de la structure qui se ruine de la manière présentée en figure 10.

Figure 10. Structure ruinée par défaut de précontrainte

 11

5.2. Coupôle de Junas

Le client Régis Deltour, membre de l’Association des Compagnons Passants
Tailleurs de Pierre, a proposé de participer à un projet de conception d’une coupôle
en pierre. Ce projet rentre dans le cadre des Rencontres de la pierre de Junas qui
rassemblent des tailleurs de pierre et des professionnels de la filière venus de France
et d’Europe.
• Contraintes et exigences
La structure est une voûte en pierre de huit mètres de diamètre et quatre mètres de
hauteur, elle s’inscrit dans une demi-sphère. Elle comporte un trou au sommet d’un
mètre de diamètre. Elle repose sur cinq piliers de deux mètres de hauteur, soixante-
dix centimètres de largeur et quarante centimètres d’épaisseur. Les blocs sont rigides
et l’interface est de type contact frottant de coefficient 0.8. Une ceinture métallique
est introduite au niveau des sommiers (modélisée par des ressorts entre les blocs)
simulant une pré-contrainte.
Pour une raideur k du ressort de l’ordre de 2OMN/m, une pré-déformation de 0,05 et
une force de pré-contrainte équivalente de l’ordre de 400kN, on obtient après 8000
itérations (pas de temps de 10-4s) :

Figure 11. Représentation de la répartition Figure 12. Représentation de la répartition

des efforts des vitesses dans les blocs

La répartition des efforts est homogène et symétrique au niveau des pieds de la

coupôle. La vitesse des blocs est de l’ordre de 10-7m/s donc négligeable : la structure
reste stable.( pour plus de détail, on pourra se reporter à R.Péralès & al, 2005)
5.3. Pont Julien
La structure est constituée de 3267 blocs. La modélisation du pont Julien reprend
la géométrie des blocs de la structure réelle (Dubois & al 2004, Chetouane 2004)
 12

L’interaction entre les blocs est de type frottement sec de Coulomb, avec une valeur
du coefficient de 0.8. La structure est modélisée sous poids propre (aucune pré-
contrainte).
Après 10000 pas de temps (pas de temps 10-4s), la structure reste stable et la
répartition des efforts est montrée ci-dessous :

Figure 13. Représentation de la répartition des efforts.

Les blocs représentés en jaune et vert sont ceux qui sont le plus chargés (les plus
contraints).

Conclusions :
Les travaux actuellement en cours concernent l’optimisation des algorithmes de
détection et traitement de contacts face-face des blocs rigides en trois dimensions. Ils
permettent d’envisager le traitement de structures réelles sous sollicitations
dynamiques, dans un proche avenir.

6. Références :

Abadie M. (1998), “Simulation dynamique de mécanismes. Prise en compte du

contact frottant”, Thèse de l’Université de Montpellier II.
Alart P., Curnier A. (1991), “mixed formulation for frictional contact problems
prone to Newton like solution methods”, Comput. Methods Appl. Mech. Eng. 92,
pp.353-375.
Baudoin A., Follenius N., Gregorowicz S. (2006), “Rapport d’étude : Escalier
Ridolfi”, Projet long de l’école des mines d’Alès.
Brogliato B., ten Dam A.A., Paoli L., Génot F., Abadie M. (2002), “Numerical
simulation of finite dimensional multibody nonsmooth mechanicals systems”,
Applied Mechanics Reviews 55(2), pp 107-150.
Cambou B., Jean M. (2001), “Micromécanique des matériaux granulaires”,
Hermes Sciences.
 13

Chetouane B. (2004), “Approche combinée éléments finis/éléments discrets pour la

modélisation des structures maçonnées”, Thèse de l’Université de Montpellier II.
Cundall P.A., Strack O.D.L. (1979), “A discrete numerical model for granular
assemblies”, Géotechnique 29(1),pp 47-65.
Dubois D., Annarumma C., Desrez E., Louis D. (2002), “Rapport d’étude :
Stabilité de structures discontinues en génie civil : application au Pont Julien”, Projet
long de l’école des mines d’Alès.
Jean M., Moreau J.J. (1992), “Unilaterality and dry friction in the dynamics of
rigid body collections”, Proc. Contact Mechanics Int. Symp., Edt A. Curnier, pp 31-
48.
Jean M. (1993), “Simulation numérique des problèmes de contact avec frottement”,
Journées de la Société Tribologique de France, Matériaux et Techniques 1-3.
Jean M. (1999), “The non-smooth contact dynamics method » Computer methods in
applied mechanics and engineering”. Comput. methods appl. mech. eng., vol. 177 ,
no 3-4 , pp. 235 – 257.
Moreau J.J. (1988), “Unilateral contact and dry friction in finite d freedom
dynamics”, CISM Courses and Lectures, 302 Springer-Verlag, pp 1-82.
Moreau J. J. (1994), « Some numerical methods in multibody dynamics :
application to granular materials » , Eur. J. Mech. A/solids, 13(4-suppl.), 93-114.
Moreau J. J. (1997), “ Numerical investigation of shear zones in granular
materials”, Proc. HLRZ-Workshop on friction, Arching, Contact Dynamics, World
Scientific, Singapore, pp233-247, (Grassberger, P. et Wolf, D. eds).
Moreau J. J. (2000) ,« Contact et frottement en dynamique des systèmes de corps
rigides ». Revue Européenne des Eléments Finis, 9, 9-28. (in French).
Péralès R., Vinches M., Bohatier C. (2005), “ERP en pierre massive. Apport des
modélisations discontinues”. VIIème rencontre de l’APS.