VuHieuNGUYEN These
VuHieuNGUYEN These
VuHieuNGUYEN These
DOCTEUR
DE
LECOLE NATIONALE DES PONTS ET CHAUSSEES
presentee par
Vu Hieu NGUYEN
Sujet de la these :
Jexprime toute ma profonde reconnaissance a Denis Duhamel, mon directeur de these, qui ma
temoigne de sa confiance et de son aide scientifique et qui, par son experience et sa competence,
ma transmis, pas a pas, sa passion pour la Dynamique des Structures. Il ma laisse une grande
liberte dans la recherche mais a ete toujours disponible pour mes questions. Sans lui, cette these
naurait surement jamais vu le jour.
Je tiens a remercier Alain Ehrlacher, le directeur du LAMI (Laboratoire Analyse des Materiaux
et Identification), pour mavoir accueilli dans son laboratoire. Jai pu travailler dans un environ-
nement de recherche exceptionnel autant du point de vue des qualites humaines que celles scienti-
fiques. Quil en soit remercie.
Je remercie vivement Pierre-Etienne Gautier, qui ma fait lhonneur de presider mon jury de
these, Didier Clouteau et Donatien Le Houedec, qui ont accepte den etre les rapporteurs et que
je remercie pour le temps quils ont passe a evaluer mon manuscrit.
Ma gratitude va a Boumedienne Nedjar pour sa gentillesse et son interet scientifique. Il a su
me consacrer du temps chaque fois que je le sollicitais.
Mes plus vifs remerciement vont egalement aux examinateurs, Gerard Gary et Karam Sab qui,
par de nombreuses discussions, par les idees precieuses et importantes, ont contribue a lorientation
de ce travail.
Je suis aussi reconnaissant au LMS (Laboratoire de Mecanique des Solides, Ecole Polytech-
nique), pour les moyens qui ont ete mis a ma disposition pour la partie experimentale.
Jadresse de chaleureux remerciements a tous les membres du LAMI, permanents, thesards ou
stagiaires, avec qui les echanges scientifiques ou amicaux mont fait de bons souvenirs inoubliables.
Je tiens a remercier du fond du coeur ma famille pour son soutien et son encouragement sans
faille malgre la grande distance. Cest son amour qui ma aide a passer les moments les plus difficiles
pendant ces annees detudes a letranger. Je remercie aussi tous mes amis pour nos discussions plus
ou moins serieuses et avec qui jai passe dagreables moments. Merci a Mai Trang pour sa patience,
sa comprehension et son encouragement.
Table des matieres
Introduction 6
1 Synthese bibliographique 11
1.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2 Les deux applications. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1 Voies ferroviaires. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1.1 Description de la structure des voies ferrees. . . . . . . . . . . . . 11
1.2.1.2 Modelisation des voies ferrees. . . . . . . . . . . . . . . . . . . . . 14
1.2.2 Voies routieres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.2.1 Modele de dimensionnement des chaussees . . . . . . . . . . . . . 16
1.2.2.2 Methode de calcul des chaussees - ALIZE . . . . . . . . . . . . . . 17
1.2.3 Sur la vibration induite par le trafic. . . . . . . . . . . . . . . . . . . . . . . 17
1.3 Rappel sur les problemes de dynamique des structures . . . . . . . . . . . . . . . . 17
1.3.1 Probleme dynamique general dans un milieu continu . . . . . . . . . . . . . 17
1.3.2 Probleme elastodynamique pour le milieu infini . . . . . . . . . . . . . . . . 18
1.3.3 Probleme de charge mobile . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4 Methodes de resolution et resultats existants . . . . . . . . . . . . . . . . . . . . . 23
1.4.1 Methode semi-analytique . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.4.2 Methode des elements finis - Condition aux limites absorbantes . . . . . . . 27
1.5 Sur le comportement dynamique du materiau ballast. . . . . . . . . . . . . . . . . 30
1.5.1 Modele discret. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.5.2 Modele continu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Annexes 159
Durant les annees recentes, les transports ferroviaires et routiers, aussi bien pour les voyageurs
que pour les marchandises, se sont beaucoup developpes avec comme consequence un trafic qui
devient beaucoup plus rapide et plus sur. Beaucoup de recherches ont ete realisees afin dameliorer
la vitesse des trains avec des voies qui sont dimensionnees pour pouvoir supporter des vitesses de
vehicules qui sont de plus en plus elevees, surtout dans les pays developpes : le TGV en France,
le Shinkansen au Japon, lICE en Allemand, lX2000 en Suede, le Pendolino en Italie, le Thalys
et lEuroStar, etc. En France, le reseau de TGV de la SNCF permet des trains qui peuvent aller
jusqua la vitesse de 300 km/h en regime commercial et de 500 km/h en essai.
Cependant, plus la vitesse des vehicules est elevee, plus le mouvement dynamique du rail et
de linfrastructure est important. La meme observation peut etre faite pour les voies routieres.
Les etudes theoriques et experimentales montrent que la vibration induite par le trafic routier ou
ferroviaire provoque des effets considerables sur le vehicule lui-meme, sur la voie et linfrastructure,
et aussi, sur lenvironnement dans la zone proche. La dynamique de structures soumises au mou-
vement des vehicules peuvent engendrer differents types de problemes : le probleme dinteraction
dynamique vehicule-structure, le probleme dynamique des voies et de linfrastructure en interac-
tion avec le sol, le probleme de propagation dondes du au passage des vehicules et ses influences
sur les batiments, etc.
La motivation initiale de cette etude est de comprendre les phenomenes physiques de vibration
apparaissant dans les voies en fonction de la vitesse des vehicules. Les mesures in situ montrent
que, lorsque le vehicule circule avec une vitesse assez elevee, surtout sur les zones ou le sol est
souple, des dommages importants peuvent etre engendres dans la voie. La deformation de meme
que le tassement de la voie sont beaucoup plus importants en tenant compte de leffet dynamique.
Plus gravement pour les voies ferrees, sous certaines vitesses de passage des trains, les ondes se
propageant dans la couche de ballast peuvent conduire a la ruine complete de la voie. Evidemment,
dans ces cas, le comportement des structures ne reste plus dans la zone lineaire mais doit etre non-
lineaire. Ces non-linearites, reciproquement, influencent le mecanisme de propagation dondes dans
la structure et ne peuvent pas etre negligees. Ceci est le deuxieme element qui motive cette etude.
Dans le cadre de cette these, le probleme se ramene a considerer le comportement dynamique
non-lineaire dune structure de forme infinie dans une certaine direction soumise a des charges mo-
biles se deplacant a une vitesse constante (ce qui peut representer des voies routieres ou ferroviaires
soumises a des passages de vehicules)
Synthese bibliographique
1.1 Introduction.
Ce chapitre presente une synthese bibliographique du probleme dynamique des structures infi-
nies soumises a des charges mobiles. On va dabord presenter une introduction generale sur les deux
applications auxquelles on sinteresse : les voies ferrees et les routes. On decrira leurs constitutions
et leurs roles, ainsi que les modeles utilises pour le dimensionnement. Puis, on va rappeler quelques
notions importantes sur la dynamique des structures (le systeme dequations dequilibre dynamique,
la definition des types dondes ...) et sur les charges mobiles qui nous serviront apres. Ensuite, les
methodes de resolution, (semi)-analytique et numerique, et les resultats existants sur le probleme
dynamique des charges mobiles seront presentes. Finalement, une synthese complementaire sur la
modelisation du comportement dynamique des milieux granulaires sera aussi donnee.
Efforts exerces sur la voie. La voie supporte une charge statique due au poids propre applique
sur lessieu et une surcharge dynamique due a linteraction dynamique entre le rail et le vehicule.
La charge statique sur un essieu est indiquee en fonction de la categorie de voie qui est definie
par lU.I.C. 1 : A : 16 tonnes par essieu, B : 18 tonnes par essieu, C : 20 tonnes par essieu.
En interaction avec les vehicules, la voie supporte des surcharges verticales, transversales et
longitudinales. Les efforts longitudinaux dus essentiellement aux accelerations et aux freinages sont
peu importants (ils peuvent cependant poser eventuellement des problemes sur certains ouvrages
dart). Les efforts verticaux et lateraux sont plus importants et causent des effets differents. La
surcharge verticale est due a limperfection de la surface du rail et aussi a la vibration de la structure
de vehicule en passant sur les traverses. Elle depend fortement de la vitesse du train (elle peut
atteindre 0.5 fois la charge statique). Les mesures montrent que les surcharges dynamiques sont
des signaux aleatoires dont la dispersion crot avec la vitesse et qui dependent de letat du rail et
aussi du vehicule. Par ailleurs, en analysant des spectres daccelerations, trois bandes de frequences
peuvent etre distinguees comme indique dans [6] :
De 0 Hz a 20 Hz : cette bande de frequence correspond a loscillation des masses suspendues du
vehicule.
1
U.I.C. : Union Internationale de Chemin de Fer.
12 1. Synthese bibliographique
blochet
190 (bton arm) entretoise
mtalique
223
72
75
260
840 735 840
172
16.5
290
445 395
1525
150
2415
ballast
4
couches o
dassise
couche "sousballast"
souscouche couche de "fondation"
nte"
couche "anticontamina
gotextile
couche de forme
plateforme
dispositif longitudinal
dassainissement
(foss ou collecteur drainant)
Figure 1.2. Profil en travers schematique des couches dassise dune voie en courbe mises
en oeuvre sur ligne nouvelle dans le cas dune plate-forme argileuse (la sous-
couche se presente alors sous forme dun multicouche)(dapres Alias [6]).
14 1. Synthese bibliographique
La recherche en mecanique des voies ferrees est toujours realisee avec deux approches paralleles :
lexperience et la simulation. Pour la simulation numerique, la premiere tache est de modeliser
la structure et aussi lexcitation. Selon le probleme auquel on sinteresse, les composantes seront
modelisees de facon plus ou moins detaillees, couplees avec les autres ou considerees independantes.
On va presenter brievement ci-dessous les modeles qui sont souvent utilises dans le dimensionnement
des voies ferrees. Un resume de ces possibilites est montre dans la figure 1.3.
couche 1
SOL couche 2
couche 3 ATTACHES
- ressort - amortisseur
- solide viscoelastique
- modele de Winkler BALLAST
- milieu multicouche infini
- ressorts simples
- modele - milieu continu
- milieu discret
- lineaire
- comportement - elastique nonlineaire
- elastoplastique
Modelisation des vehicules. Dans la litterature, la structure des vehicules est souvent modelisee
par un systeme multi-essieu complexe de masses-ressorts-amortisseurs. Par exemple, on peut trou-
ver dans [46] la structure dune voiture de train (qui est constituee par la caisse, la suspension
primaire et secondaire, le boggie, les essieux). Dans ce modele, la caisse et le boggie sont consideres
comme des corps rigides. Ce modele est utilise pour determiner la surcharge dynamique sur la
voie, voir [6] [46]. Dans une etude plus detaillee, on utilise des calculs dynamiques du probleme
dinteraction entre le vehicule et la voie [63] dans lequel la surface irreguliere du rail peut etre
introduite [39][96].
Le calcul des structures de voies peut se simplifier en considerant le vehicule comme une charge
ou un groupe de charges. Les hypotheses sur les charges sont variees et dependent du niveau de
complexite du probleme pose : elles peuvent etre fixes ou mobiles, lamplitude peut etre constante,
harmonique, variable en fonction du temps, aleatoire, etc.
Modelisation des structures. On ne presente ici que les modeles qui concernent la modelisation
globale des voies ferrees.
Les rails. Le modele le plus simple (et aussi le plus utilise) est la poutre dEuler-Bernouilli
[6][39], etc. Lhypothese de Timoshenso est aussi utilisee pour tenir compte plus precisement de
linfluence du cisaillement [47][21], etc.
Les semelles. Elles ne sont pas toujours prises en compte dans les calculs. Si on doit en tenir
compte, elles peuvent etre considerees comme un ressort-amortisseur [6] [101, 100].
Les traverses. Elles sont assez rigides par rapport aux autres composantes. Elles peuvent etre
une masse ponctuelle dans un calcul 1D ou un corps rigide de 3 ou 6 ddl selon les cas 2D ou 3D.
Dans le cas 3D, les traverses de 2 blochets sont liees par une poutre qui represente les entretoises
[16].
Le ballast. Cest un point tres difficile dans la modelisation des voies ferrees a cause de ses
proprietes discretes. Pour un modele 1D plus simple, il peut etre modelise par un systeme continu
de ressorts dont la masse est uniforme [101, 100]. Les problemes en 2D ou en 3D proposent des
modeles soit continus (elastique [12][16] ou non-lineaire [8],[37]...) soit discrets [103, 104], etc.
Une discussion plus detaillee sur le ballast sera presentee plus loin dans la section 1.5.
La plate-forme et le sol. Les modeles 1D considerent linfrastructure comme un systeme de
ressorts, dit la fondation de Winkler. Dans les calculs en 2D ou 3D, ils sont souvent modelises
par un milieu multicouche infini. Le comportement du sol peut etre lineaire ou non-lineaire en
fonction du type de materiaux (Tresca, Coulomb-Drucker ...)
Le modele global de calcul est une combinaison des differentes composantes et il existe plusieurs
facons possibles de modeliser la voie ferree, de la plus simple a la plus complexe.
Quelques etudes existantes. On peut citer ici quelques modedes existants pour le dimension-
nement des voies ferrees. Balsan [12] a realise une analyse en 3D par elements finis du comportement
de la structure dassise dans les essais de Vienn-Asernal et de Derby. Sauvage [99] et ensuite Pro-
fillidis [94] ont developpe en utilisant le code Rosalie du LCPC des modeles assez complexes (rail,
traverses, blochets, semelles) pour etudier le comportement des voies soumises a des efforts ver-
ticaux et statiques. La loi de plasticite de Coulomb-Drucker a ete utilisee pour la couche dassise
(dont les parametres de critere de plasticite sont determines par des essais triaxiaux). Un autre
modele, appele modele RB3D a ete developpe par le CERAM (Alaoui et Naciri (1995) [4], Guerin
[46]). Il permet devaluer les evolutions des defauts de voie a partir de la charge des essieux (avec
la surcharge dynamique dependante de la vitesse), de la nature du sol et de la loi de tassement
du ballast. Recemment, Bodin a etendu ce modele pour obtenir le modele RB3L qui permet de
considerer aussi linfluence des efforts lateraux [16][17].
16 1. Synthese bibliographique
Il nexiste encore que tres peu detudes dynamiques pour le dimensionnement des voies bal-
lastees. Aubry et al. [8] ont calcule en utilisant le logiciel MISS3D (qui est developpe a lEcole
Centrale de Paris) une structure de voie ou la couche de ballast est consideree comme un materiau
elasto-plastique. Le couplage des elements finis-elements de frontieres a ete utilise. Les tassements
residuels sont compares dans deux cas de comportement : elastique et elastoplatique.
Les sols en place etant incapables de supporter les charges induites par le trafic, la chaussee a
pour role de diminuer, grace a lapport de couches de materiaux, les contraintes sur le sol-support.
Une chaussee est donc constituee dun empilement de couches qui a deux composantes principales :
le corps de chaussee et la couche de roulement.
Le corps de chaussee. Il est appele aussi couches dassises et permet de repartir les charges
induites par les vehicules afin dobtenir des contraintes acceptables pour le sol-support. Il est
constitue de deux sous-couches, appelees couche de base et couche de fondation.
la couche de roulement. Son role est de proteger les assises des infiltrations deau et de lagressivite
du trafic.
Deux couches supplementaires peuvent eventuellement etre ajoutees pour remedier a des problemes
specifiques : la couche de liaison qui assure un collage entre la couche de roulement et le corps de
chaussee et la couche de forme qui reprend les heterogeneites trop importantes du sol-support [3].
d = 3a
a = 12, 5 cm
a
q0 = 0, 662 M P a
q0
Couches de surface
couche de roulement
couche de liaison (eventuelle)
couche de base
Corps de chaussee
couche de fondation
couche de forme (eventuelle)
sol
de rayon a = 12.5 cm avec une distance entre elles d = 3a. La pression est uniformement repartie
avec la valeur q0 = 0.662 M P a [13].
Au LCPC, le calcul des contraintes et des deformations dans cette structure est effectue par
ALIZE, qui est un programme de dimensionnement des chaussees. Ce programme permet de donner
les champs statiques de deplacement et de contraintes dans la structure soumise a une ou a un
groupe de charges circulaires appliquees sur la surface libre. Il est base sur le modele de Burmister
(1943) dans lequel les couches sont supposees infinies et planes ; la charge est circulaire et le
probleme se traite par consequent en coordonnees cylindriques grace a la symetrie. Brievement, la
resolution seffectue de la facon suivante : dans chaque couche, le deplacement et les contraintes
sont representes par les derivees de la fonction potentielle, ensuite, en utilisant la transformation de
Hankel, il apparat une equation differentielle dune seule variable dans la direction verticale dont
on peut presenter la solution generale avec 4 coefficients inconnus. Ces 4 inconnues seront trouvees
a partir des conditions de continuite entre les couches a la frontiere. Finalement, la transformation
inverse est calculee numeriquement par la methode de quadrature de Gauss pour obtenir la solution.
Les vibrations induites par les mouvements des vehicules routiers ou ferroviaires sont des effets
tres significatifs sur lenvironnement, le sol, les ouvrages dans la zone voisine, et aussi sur la
stabilite des vehicules et des structures eux-memes. Les etudes de vibration induites par le trafic se
divisent en deux problemes principaux : probleme demission directe de bruit (bruit de roulement)
et probleme de propagation dondes dans la structure.
Emission directe de bruit. Le premier sinteresse a la radiation du bruit emis directement a
partir de linteraction entre les roues des vehicules et le rail. Le but est de predire et de donner
des indications sur la reduction du bruit. Beaucoup de recherches dans ce domaine sont realisees
sur le logiciel TWINS (Track-Wheel Interaction Noise Software) qui est developpe par lEuropean
Rail Research Institute (ERRI) [45][55].
Propagation dondes dans la structure. Le deuxieme sinteresse aux ondes qui ne sont pas
emises directement dans lair comme le bruit mais qui se propagent dans les infrastructures et le
sol au-dessous. Le but est devaluer la reponse dynamique des vehicules - structures et de predire
les criteres de ruine dans le dimensionnement des voies ferrees ou des chaussees.
Dans le cadre de cette these, nous nous interessons au deuxieme probleme.
f :
( , t). = ( , t)
u : =
Figure 1.5.
Conditions initiales
E E
= , = (1.8)
(1 + )(1 2) 2(1 + )
Le probleme dynamique presente au dessus dans ce cas elastique peut se representer par
lequation de Navier [2, 38] comme suit :
2 ui (x, t) + ( + ) ( u(x, t)) = ui (x, t) (1.9)
xi
2 2 2
ou 2 designe loperateur Laplacien : 2 = 2 + 2 + 2 et est la notation de loperateur
x1 x2 x3
u1 u2 u3
de divergence : u = + + .
x1 x2 x3
1.3 Rappel sur les problemes de dynamique des structures 19
Ondes volumiques. En absence de condition aux limites, la solution generale des champs de
deplacement u de lequation de Navier peut etre decomposee dans les deux types dondes de natures
differentes suivantes :
ou les fonctions (x, t) (scalaire) et (x, t) (vectorielle) sont appelees les deux potentiels de Lame
et verifient les equations des ondes :
s
2 + 2
c2p 2 2 = 0, cp = (1.11)
t
r
2
c2s 2 2 = 0, cs = (1.12)
t
Cela montre deux mecanismes independants de propagation donde a deux vitesses finies :
onde de compression (appelee aussi onde de dilatation ou onde irrotationnelle) associee au
mouvement des particules qui seffectue parallelement a la direction de propagation de londe
et a la vitesse cp . Cette propagation est due aux efforts de traction-compression induits par
leloignement ou le rapprochement relatif des faces perpendiculaires au mouvement et est notee
londe P.
onde de cisaillement (appelee aussi onde de distorsion ou onde equivolumique) associee au
mouvement des particules qui seffectue perpendiculairement a la direction de propagation de
londe et a la vitesse cs . Ces ondes se propagent plus lentement, sous leffet des efforts de ci-
saillement induits par le glissement relatif des faces perpendiculaires au mouvement. Ses deux
composantes dans ce plan sont notees londe SV et londe SH.
q
12
Notons que le rapport entre cp et cs ne depend que du coefficient de Poison : c s = cp 2(1) ce
qui montre quon a toujours cs < cp car > 0.
B B
onde SV
A A onde SH
Onde de surface. La presence dune surface libre perturbe les ondes simples de depart par
lapparition de phenomenes de reflexion et de refraction. Sur la surface libre dun demi-espace
infini, il apparat un troisieme type donde dont lamplitude est maximum a la frontiere et diminue
exponentiellement avec la distance a surface libre. Celle-ci est appelee londe de surface ou londe
20 1. Synthese bibliographique
de Rayleigh 2 car la perturbation est pratiquement confinee dans une couche mince adjacente a la
surface. Dans le cas dun demi-espace elastique, londe se propage a une vitesse c R < cs qui est
determinee par la racine de lequation de Rayleigh :
1 1
R(cR ) = [2 (c2R /c2s )]2 + 4[1 (c2R /c2p )] 2 [1 (c2R /c2s )] 2 = 0 (1.13)
Le mouvement des particules pour ce type donde est de nature elliptique et retrograde par rapport
a la direction de propagation de londe.
Si la surface libre nest par parfaitement plane, il peut apparatre une autre onde de surface
(dite londe G, dapres Le Tallec [71]) qui est guidee localement par la direction de la discontinuite.
Elle se propage a une vitesse vG qui est inferieure a cR .
Les energies de ces ondes de surface se concentrent localement sur la surface et en plus,
se dispersent tres lentement (notons aussi quen absence damortissement, les ondes volumiques
sattenuent a linfini avec une decroissance de lordre de 1/r 2 mais de lordre de 1/ r pour celle
de Rayleigh, r est la distance entre la source et le point considere). Les ondes de surface sont les
causes principales des risques lors dun seisme.
direction de propagation
Ondes dans un milieu stratifie. En realite, le sol ne peut pas etre modelise par un massif
semi-infini mais plutot par un milieu stratifie (ou milieu multi-couche) qui est constitue de plu-
sieurs couches. En presence dune interface entre deux couches de proprietes physiques differentes,
une onde incidente voit son type modifie ainsi que sa direction de propagation et son amplitude
du fait des phenomenes de refraction et reflexion. Pour le cas plus simple correspondant a deux
demi-espaces en contact, quand une onde P rencontre linterface, elle donne deux ondes P et SV
reflechies, et deux autres ondes P et SV refractees ; une onde SH incidente ne donne que des ondes
SH refractees et reflechies [2, 38].
La reponse dynamique complete dun domaine elastique a plusieurs couches devient rapidement
beaucoup plus difficile a determiner. Le mouvement dans une couche est une superposition des
ondes reflechies et refractees qui viennent des interfaces de discontinuite. En outre, des ondes de
surface apparaissent aussi a chaque interface. Le calcul de la reponse dynamique se fait alors par
integration numerique ou par analyse en modes propres.
2
Le probleme de londe de surface dun demi-espace elastique a ete etudie par Rayleigh en 1885.
1.3 Rappel sur les problemes de dynamique des structures 21
Condition de radiation. Dans le cas ou le domaine considere est infini, la condition de radiation
du champ de deplacement dun type donde (dont la vitesse est c) pour la frontiere infinie secrit :
s1 u 1 u
lim r 2 + 0 , lim u 0 (1.14)
r r c t r
Notation. On sinteresse au probleme ou la charge se deplace avec une vitesse constante v sur
une ligne droite dans le domaine (ou sur sa surface libre). La charge est soit constante soit variable
en fonction du temps (harmonique, par exemple). On peut definir la forme generale dune charge
mobile sur une ligne e comme suit :
ou f 0 (t) est le vecteur damplitude qui ne depend que du temps, est la fonction de Dirac.
Repere mobile. Dans les etudes des problemes de charge mobile avec une vitesse constante,
plusieurs auteurs (Fryba [39], Eringen et Suhubi [38], etc) preferent passer dans un repere mobile
qui sattache a la position de la force. Ce fait est realise par un changement de variable. Par
exemple, si la charge se deplace suivant laxe Ox avec la vitesse v, on a le changement de variable
suivant :
x = x vt ; y = y ; z = z (1.17)
Le probleme est reecrit dans ce nouveau repere avec les variables x , y , z . Dans certains cas, par
exemple si la charge est constante ou harmonique, le probleme se simplifie beaucoup dans le repere
mobile.
repere fixe repere mobile
z z
O x O x
y y
linstant 0 linstant t
vt
Regimes dune charge mobile. En analysant le probleme de propagation dondes avec une
source mobile, il est important de considerer le rapport entre la vitesse de la source v et celles des
ondes du milieu. Le probleme ou la source est en mouvement a ete considere en premier pour les
22 1. Synthese bibliographique
phenomenes acoustiques. Dans un domaine elastique homogene, on introduit trois notions, dites
nombres de Mach, qui representent le rapport entre v et c p , cs et cR :
v v v
Mp = ; Ms = ; MR = (1.18)
cp cs cR
Pour le probleme elastodynamique dun milieu massif, on peut distinguer trois regimes differents
dependant de ces nombres de Mach [38] :
Regime subsonique. Quand Ms < 1 : la vitesse de la source est inferieure de celle de londe S.
Regime transonique. Quand Ms > 1 > Mp : la vitesse de la charge est entre celle de londe S et
celle de londe P .
Regime supersonique. Quand Mp > 1 : la vitesse de la charge est superieure celle de londe P .
Cone de Mach
cs (t ) de londe S Cone de Mach Cone de Mach
cs (t ) de londe S de londe P
cp (t ) cp (t ) cs (t )
cp (t )
A B C A B C A B C
cp t cp t cp t
cs t cs t cs t
La figure 1.9 presente une illustration de ces 3 regimes : sub-, trans- et super-sonique. Elle montre
la propagation des ondes quand on regarde une source qui part de A a C en passant par B
(AB = v < AC = vt). Quand v depasse cs ou cp , les surfaces de londe de choc, dites les cones
de Mach, apparaissent.
En outre, concernant londe de surface, il sera utile de comparer la vitesse de la charge vis a
vis de la vitesse donde de Rayleigh. En considerant la valeur du nombre de Mach M R , on peut
distinguer 2 regimes differents :
si MR < 1 : regime sub-Rayleigh.
si MR > 1 : regime super-Rayleigh.
2 2 2 2v 1
(1 Mp2 ) + 2+ 2+ 2 = 0 (1.19)
x21 x2 x3 cp x1 c2p
2 2 2 2v 1
(1 Ms2 ) 2 + 2 + 2 + 2 2 = 0 (1.20)
x1 x2 x3 cs x1 cs
1.4 Methodes de resolution et resultats existants 23
Il est souligne dans [38] que la nature de ces equations en fonction des variables despace change
suivant les valeurs des nombres de Mach M : quand M < 1, on a une equation elliptique et quand
M > 1, on a une equation hyperbolique.
est donnee dans [38] a partir des solutions du probleme de Lamb 3 . Les calculs montrent que la
premiere resonance apparat quand la vitesse de la charge sapproche de la vitesse de londe de
Rayleigh.
A partir de ces demarches, des methodes de calcul semi-analytiques ont ete developpees dans les-
quelles les transformations inverses sont realisees numeriquement. Celles-ci permettent de resoudre
des problemes plus compliques, quand on ne peut pas avoir une forme analytique explicite des
solutions dans le domaine frequentiel, par exemple quand le domaine nest plus homogene mais est
multicouche. Pour evaluer la transformation de Fourier inverse, la technique de transformation de
Fourier rapide (FFT) est beaucoup utilisee.
Dans [27][28], Barros et Luco ont presente la methode pour obtenir les deplacements et les
contraintes en regime permanent dans un milieu multi-couche visco-elastique. La solution est ob-
tenue par la transformation de Fourier inverse qui est calculee en utilisant la technique de FFT et
lintegration numerique adaptative de Filon. La solution est validee en comparant avec la solution
analytique dans le cas 2D.
Une etude detaillee du modele bidimensionnel (homogene ou multi-couche) est presentee dans
la these de Leufeuve-Mesgouez [73]. Ladimensionnement des resultats par rapport a la longueur
donde de Rayleigh a ete presente. En plus, un calcul analytique approche qui prend en compte
uniquement londe de Rayleigh a ete aussi propose.
Le probleme dun demi-espace homogene soumis a une charge uniforme mobile distribuee sur un
rectangulaire peut etre trouve dans les travaux de Jones et al. [60] (pour une charge constante), [58]
(pour une charge harmonique). Leffet de plusieurs charges a ete discute dans [72] ou les resulats
de deplacement dus a un groupe de 9 charges rectangulaires posees a des distances differentes
montrent lexistance dinterferences constructives ou destructives. Picoux et Le Houedec [92] ont
etudie leffet dun effort lateral et vertical sur une charge harmonique et ont montre que cet effet
nest pas negligeable pour les deplacements lateraux.
Dans le cas dun milieu multi-couche soumis a une force constante, la vitesse critique, qui etait
egale a la vitesse donde de Rayleigh (qui ne depend que des parametres physiques (E, )) depend
aussi des parametres geometriques, i.e des epaisseurs des couches. On peut trouver dans [59], [73] les
etudes sur la couche critique du probleme dune couche elastique posee sur un demi-espace rigide
ou elastique. La vitesse critique due a une charge harmonique a ete etudiee par Dietermann et
Metrikine [31]. Lepaisseur critique de la couche elastique est aussi determinee dans cette reference
en fonction de la vitesse dans le cas ou la frequence de la charge est proportionnelle a sa vitesse
(quand la vitesse du train est grande, la frequence de la charge est notamment determinee comme
un rapport de la vitesse et de la distance entre les traverses). En general, quand la profondeur de
la couche elastique augmente, la vitesse critique est diminue.
Quand la vitesse de charge depasse la vitesse critique, la solution dans le domaine frequentiel-
nombres dondes devient tres singuliere et celle-ci pose des difficultes pour evaluer les integrations
de la transformation inverse. Du point de vue analytique, on peut faire lintegration dans le plan
complexe ou on doit specifier les poles, choisir les contours et calculer les residus [38][39][73] ... En
revanche, les points singuliers ne doivent pas apparatre dans une solution numerique. Plusieurs
auteurs ont introduit une petite viscosite pour deplacer des points singuliers de laxe dintegration
afin davoir une solution numerique finie. Neanmoins, la fonction a integrer est tres irreguliere dans
ce cas et un nombre de points de calcul assez important est demande pour evaluer precisement
lintegration [58]. La technique dintegration numerique adaptative a ete aussi utilisee [27][54]. Par
ailleurs, la methode de transformation par ondelettes qui permet dutiliser un maillage dintegration
adaptatif en localisant les points singuliers est proposee par Grundmann et al. [44], Lieb et Sudret
[75] ...
3
Lamb (1904) a resolu le probleme transitoire dun demi-espace soumis a une force ponctuelle ou a une ligne de
force a la surface libre : f (t) = f0 (x)(t).
1.4 Methodes de resolution et resultats existants 25
Dans les etudes quon a montrees, la transformation de Fourier est tres souvent utilisee car la
plupart etaient interessees par le regime permanent du probleme. Si on sinteresse a la solution
transitoire, la transformation de Laplace peut etre utilisee. Cette methode a ete proposee par
Gakenheinmer et Miklowitz (1969) et completee par Bakker et al. [11] afin detudier le probleme
dun demi- espace soumis a une force ponctuelle mobile.
Poutre infinie - fondation de Winkler. Ce systeme est le plus simple pour modeliser le
rail couple avec la fondation. Il etait premierement etudie par Timoshenko depuis lannee 1926.
La methode de resolution est toujours similaire a celle du cas dun massif infini sauf quon doit
resoudre une equation differentielle du quatrieme ordre. Dans [39], Fryba presente en detail un
calcul analytique du probleme dune force constante se deplacant sur une poutre dEuler-Bernouilli
infinie sur la fondation de Winkler pour toutes les possibilites de vitesses et de valeurs de viscosite.
Une vitesse critique a laquelle la resonance apparat a ete trouvee, elle depend de la densite et de la
rigidite en flexion de la poutre et aussi de la rigidite de la fondation. La vitesse critique augmente
quand la rigidite de la fondation augmente.
Lhypothese dEuler-Bernoulli peut ne pas etre suffisante dans certains cas. Chen et Huang [21]
ont presente un calcul semi-analytique pour tous les cas de vitesse pour une poutre de Timoshenko
soumise a une force ponctuelle constante. Ils ont montre que quand la rigidite de la fondation
devient importante (i.e leffet de cisaillement dans la poutre est plus important), la vitesse critique
calcule avec lhypothese dEuler-Bernoulli est superieure a celle de Timoshenko.
Si la charge est harmonique, les ondes qui sont dues a la vitesse et aussi a la frequence de
la charge peuvent se propager de differentes manieres. Dans [7], Andersen et al. ont presente
une etude adimensionnelle pour des mecanismes differents de propagation dondes dependant des
valeurs relatives de vitesse et de frequence de charge vis a vis de celles critiques.
Hardy [48] a presente une methode graphique pour identifier la phase et la vitesse de groupe
donde en utilisant une integration de convolution mobile. Lexistence et le nombre des ondes
propagatives dependent fortement de la vitesse et de la frequence de la charge.
La vibration dune masse ponctuelle se deplacant sur une poutre - fondation de Winkler en
compression a ete consideree dans [80]. Les auteurs montrent que pour une vitesse suffisamment
grande de la masse, une instabilite de vibration apparait dans la poutre. La vitesse critique est
plus petite avec une plus grande masse ou avec une plus grande force axiale.
Couplage de poutre - massif infini. Dans une etude plus detaillee, le systeme de ressorts de
Winkler ne peut plus representer la fondation reelle. Il sagit de poser un probleme de couplage
entre une poutre et un massif elastique en presence dune charge mobile. Dieterman et Metrikine
[33] ont calcule la rigidite equivalente dun demi-espace en interaction avec une poutre dEuler-
Bernoulli. Les auteurs ont montre que la rigidite equivalente depend fortement de la frequence et
du nombre donde de la poutre si la vitesse de phase est au voisinage de la vitesse de Rayleigh.
Deux vitesses criques sont aussi determinees dans le cas dune charge constante : une egale a la
vitesse donde de Rayleigh et lautre qui est un peu plus petite. Quand la vitesse de la charge est
egale a la vitesse critique, la rigidite equivalente est zero, cest a dire quil ny a pas dondes qui
sont generees dans le milieu massif.
Les memes auteurs ont presente dans [32] la reponse dynamique de ce probleme avec une charge
constante dans tous les cas de vitesses. Dans [83], la reponse dune poutre en compression posee
sur un demi-espace due a une charge laterale a ete aussi consideree. La vitesse critique dans ce cas
est determinee en fonction de la force de compression axiale dans la poutre et on peut montrer
quelle est toujours inferieure a celle avec la charge verticale.
26 1. Synthese bibliographique
Un modele plus complexe dune voie ferree constituee dune poutre (pour le rail), dun systeme
masse-ressorts-amortisseurs (pour les blochets, les semelles et le ballast) et dun demi-espace multi-
couche a ete etudie avec une charge constante ou harmonique par Sheng et al. [101][100]. Les effets
de vitesse et de frequence de charge, ainsi que du type dinfrastructure de voie sont presentes dans
le travail de Jones et al. [57].
En comparant la reponse dynamique dune charge mobile appliquee sur une poutre posee sur
une fondation de Winkler et posee sur un demi-espace elastique 3D, Dinkel et al. [34] ont propose
une methode pour determiner des parametres de Winkler qui peuvent representer le domaine 3D.
Une procedure doptimisation numerique est aussi proposee pour minimiser lerreur.
Dans [82], Metrikine et Vrouwenvelder ont considere la vibration a la surface dun milieu 2D
generee par une charge mobile sur une poutre dans un tunnel. Le modele 2D est utilise avec trois
types differents de charge : constante, harmonique et aleatoire stationnaire. Il est montre que la
vitesse minimale de phase dans la poutre est inferieure a la vitesse donde de Rayleigh dans la couche
et est plus grande avec une profondeur plus grande. Pour une charge deterministe, le deplacement
a un point dobservation a la surface est analyse par le spectre damplitude de vibration a ce point.
Pour une charge aleatoire, cest la variance de la vibration qui est consideree.
Knothe et Wu [62] ont etudie le comportement dynamique dune voie posee sur un demi-espace
ou sur une fondation viscoelastique simple. Les traverses discrets sont tenues en compte. Ils ont
observe que pour les frequences inferieures a 250 Hz, il y a des grandes differences entre les resultats
du modele dun demi-espace et du modele de fondation visco-elastique. Au-dela de cette valeur, la
difference est negligeable, le modele simple visco-elastique peut etre utilise.
Poutres sur des appuis periodiques. Dans [14] [15], Belotserkovskiy a etudie la vibration
due a une charge mobile harmonique posee sur des appuis elastiques periodiques. Des conditions
speciales sont imposees sur deux points quelconques qui sont separes par une distance qui est egale
a celle entre deux appuis. Ces conditions imposent que, dans le regime permanent, il est suffisant
de considerer uniquement le probleme dans une partie qui se trouve entre deux appuis.
En developpant le probleme presente dans [33], Metrikine et Popp [81] ont determine la rigidite
des ressorts equivalents qui peut etre mis a chaque appui et qui permet de decrire exactement
la reaction en regime permanent du demi-espace elastique. Cette rigidite (qui est une variable
complexe) est calculee en fonction de la frequence de vibration de la poutre et du dephasage de
vibration de lappui voisin. Le dephasage depend de la frequence et de la vitesse de la charge, et de
la periode dappuis de la poutre. La vitesse critique (quand la rigidite est egale a zero) est evaluee
par des relations entre le dephasage et la frequence de la poutre.
Krzyzynski et al. [70] ont etudie un probleme dans lequel le rail est modelise par une poutre
de Timoshenko et les blochets sont consideres comme des corps rigides de 2 degres de liberte. Le
theoreme de Floquet est utilise.
Pour les structures dont les geometries sont plus complexes ou avec des materiaux non-lineaires,
il faut utiliser les methodes numeriques comme la methode des elements finis. En general, lanalyse
par elements finis de la propagation des ondes dans un domaine infini pose deux problemes majeurs.
Premierement, la structure doit etre discretisee de facon telle que les ondes peuvent se propager.
Deuxiemement, comme il est seulement possible de modeliser la structure dans une partie finie, des
frontieres artificielles doivent etre introduites pour permettre de decrire linfluence des domaines
exterieurs.
Le calcul des structures soumises aux charges mobiles par la methode des elements finis a attire
beaucoup de travaux depuis une dizaine dannees. Il y a deux points de vue qui sont le plus souvent
utilises :
(i) calcul dans le repere fixe avec une source mobile
(ii) calcul dans le repere mobile avec une source fixe.
Calcul dans le repere fixe. Habituellement, lequation dynamique (et eventuellement non-
lineaire) en elements finis est etablie en discretisant la formulation variationnelle de Galerkin pour le
systeme presente dans le paragraphe 1.3.1. Ce fait nous donne un systeme dequations differentielles
en temps avec des conditions initiales :
la difference est que la position de la force dans ce cas change a chaque pas de temps. Les methodes
classiques comme la methode de superposition modale ou les procedures dintegration numerique
(-Newmark, Wilson-, -HHT ...) peuvent etre utilisees. La premiere nest valable que pour les
cas elastiques. Si la structure est non-lineaire, il faut utiliser la deuxieme pour trouver la solution
en regime permanent apres le calcul transitoire. La difficulte de cette facon de calculer vient de
deux raisons essentielles. Dabord, comme le regime stationnaire ne peut etre atteint quapres un
certain temps de calcul, i.e. apres une certaine distance de deplacement de la charge, la taille du
maillage necessaire doit etre grande. Puis, la densite de maillage au voisinage de la charge (ou la
non-linearite peut apparatre par exemple) doit se deplacer avec la charge, et donc, on a besoin
dune grande zone avec des elements fins. En plus, les schemas dintegration numerique demandent
des pas de temps qui sont suffisamment petits pour assurer la stabilite.
En utilisant les elements finis lineaires dans le calcul dynamique, on a besoin au moins de 4
elements sur une longueur donde afin de pouvoir representer le champ dondes. Cependant, pour
une bonne description des ondes, le nombre delements necessaires est de lordre de 10 (voir par
exemple [51][115]). Si des elements avec une fonction dinterpolation dordre superieur sont utilises,
le nombre delements peut etre reduit.
Dans ses etudes du probleme de poutres finies soumises a une charge mobile, Rieker et al. [97]
ont montre quon a besoin dun nombre delements pour le probleme de charge mobile qui est de
2 a 8 fois plus que dans un probleme statique. Le vecteur force equivalent pour un element qui
a une charge mobile a linterieur depend beaucoup de la fonction dinterpolation. Les problemes
1D finis (structures de poutres simples ou multi-travees couplees avec un systeme qui modelise
les vehicules) peuvent se trouver dans plusieurs etudes differentes. On peut citer par exemple les
travaux de Delgado [30], de Henchi et al. [49], de Olsson [87], de Wu et al. [109], de Zibdeh et al.
[114], etc.
Dans les annees recentes, avec lavancee des capacites des ordinateurs qui permettent de calculer
des modeles avec de gros maillages, les problemes en 3D peuvent etre abordes. Dans ces modeles,
28 1. Synthese bibliographique
la structure detaillee de la voie : rail, blochet, semelles et aussi de la structure dassises : ballast,
sol avec la non-linearite peuvent etre pris en compte.
Dans [36], Ekevid et al. ont realise des calculs transitoires en 2D non-lineaire par elements finis
en utilisant une technique de temps-espace adaptative qui est basee sur la procedure de Galerkin
discontinue (Discontinuous Galerkin 4 ). Les pas de temps et le maillage sont redefinis a chaque
instant de calcul a partir dune estimation derreur. Le cas ou on est en regime supersonique a ete
aussi traite. Un autre calcul 3D de la structure de voie ferree en non-lineaire (de comportement
elasto-plastique) est presente par les memes auteurs dans [37]. Les elements finis couples avec des
elements aux frontieres (Scaled Finite Element) sont utilises. Cependant, ces calculs demandent
des ordinateurs tres puissants et des temps de calcul enormes.
On peut aussi trouver dans [90] un calcul dynamique transitoire dun couplage nonlineaire de
chaussee - sol propose par Pan et al. . Une modelisation mixte des elements finis (pour la chaussee)
- elements frontieres (pour le sol) a ete utilisee. Un concept interessant appele pass-over a ete
propose par les auteurs pour sadapter au probleme de charge mobile : le champ de deplacement
et de contrainte se deplace avec la position de la charge apres chaque pas de temps. Les solutions
obtenues dans le pas de temps precedent sont utilisees comme les valeurs initiales pour le pas de
temps actuel.
Par ailleurs, en supposant que le domaine considere est homogene dans la direction du mou-
vement de la charge, un modele appele modele 2.5D a ete propose. Dans ce modele, une transfor-
mation est realisee dans la direction du mouvement et le domaine du plan dans les deux autres
directions est discretise par elements finis. On peut citer par exemple le travail de Yang et Hung
[110] ou les auteurs ont valide ce modele dans tous les cas subsonique, transonique et supersonique.
Un modele plus complique de couplage lineaire entre des vehicules, les infra-structures et le
sol infini a ete utilise par Clouteau et al. [23] en utilisant la technique de sous-structuration et
des elements finis aux frontieres. Les hypotheses supposent que les infrastructures sont invariantes
dans la direction de mouvement Ox des vehicules qui ont des vitesses constantes et que le sol est
stratifie horizontalement. Le probleme a ete dabord transforme par la methode integrale le long
de la direction infinie : transformation en t et en Ox dans linfrastructure et transformation de
Hankel dans le sol. Ensuite, il est resolu par la methode de Galerkin dans le domaine borne. Un
autre modele numerique qui est construit de la meme facon peut etre trouve dans [76] dans lequel
la solution numerique a ete validee en comparant avec celle analytique dans le cas dune charge
ponctuelle se deplacant sur un demi-espace homogene. Les auteurs montrent aussi les resultats
numeriques dun probleme de passage dun vehicule (qui est modelise comme un systeme 2D de 4
ddl) sur une chaussee en interaction avec le sol.
Calcul dans le repere mobile. Dans le but de diminuer la taille du domaine de calcul qui doit
etre tres grande pour une charge mobile, lidee decrire la formulation des elements finis dans un
repere mobile a ete proposee. Le probleme est etabli dans ce nouveau domaine par un changement
de variable. Lavantage de cette methode est la facilite pour determiner la reponse de la structure
en regime permanent. Cette methode est encore tres peu utilisee dans la litterature.
Andersen et al. [7] ont presente un calcul par elements finis dans le repere mobile du probleme
dune poutre dEuler-Bernoulli posee sur la fondation de Winkler. Le systeme dequations dyna-
miques de meme forme que lequation 1.21 est retrouve mais avec des termes convectifs ajoutes
dans les matrices de rigidite et damortissement qui viennent du changement de variable. Une
charge damplitude variable en fonction du temps a ete consideree. Les resultats numeriques et
analytiques sont compares dans le cas de charges harmoniques.
4
la procedure DG discretise par des elements finis le temps et lespace en meme temps avec des fonctions de forme
qui sont continues en deplacement mais peuvent etre discontinues au niveau du temps [52][53][74]
1.4 Methodes de resolution et resultats existants 29
Si la charge consideree est constante, le probleme dans le repere mobile en regime permanent
devient statique. On peut trouver dans [64] une etude sur le probleme 2D ou dans [85] une autre
sur un demi- espace non-lineaire en 3D.
Quand la vitesse de charge devient grande, les termes de convection deviennent importants et la
discretisation de Galerkin standard donne des instabilites de vibration. Dans ce cas, la methode de
Petrov-Garlerkin (voir par exemple [115]) qui introduit une partie asymetrique dans la fonction de
forme peut etre utilisee. Dans [65], Krenk et al. propose une methode de Taylor-Galerkin modifiee
ou un terme de correction de type upstream difference est implemente dans la discretisation
standard.
Cependant, la resolution par elements finis dans le repere mobile est toujours limitee pour
les vitesses subsoniques. Dans le cas supersonique, Kok [64] a propose que, au lieu decrire deux
conditions aux limites, on doit ecrire une condition initiale a la frontiere droite. De toute maniere,
cette methode ne peut sappliquer quaux cas simples. Il ny a pas encore de procedure generale
pour toutes les vitesses de charge.
Condition aux limites. En utilisant la methode des elements finis, le maillage est limite a un
domaine fini avec des conditions imposees aux limites, soit pour les deplacements soit pour les
contraintes. Cela pose toujours des difficultes pour modeliser le probleme de propagation dondes
dans un domaine infini a cause des ondes reflechies aux frontieres. En regardant dans la litterature,
on peut trouver un grand nombre detudes pour resoudre cette difficulte. Le plus important est
dassurer la condition de radiation qui impose quil ny a pas denergie venant de linfini. On peut
enumerer quelques modeles differents appliques aux problemes de propagation dondes :
Calculer sur un domaine tres grand tel que toutes les energies peuvent etre amorties avant
quelles arrivent aux frontieres. Cette modelisation est tres chere et ne peut sappliquer quaux
problemes ou les ondes sont attenuees rapidement par les amortissements dans la structure.
Utiliser des ressorts et des amortisseurs (visco-elastiques) aux frontieres tels que les frequences
dominantes sont amorties.
Utiliser des elements infinis speciaux qui peuvent sadapter aux cas specifiques. En principe, les
elements infinis sont construits par des fonctions dinterpolation speciales qui considerent les
noeuds exterieurs comme etant a linfini [115]. Pour le probleme de propagation dondes, ces
fonctions dinterpolation doivent comporter aussi des termes qui representent les ondes qui se
propagent vers linfini (voir par exemple Yerli et al. [111]).
Les methodes mixtes qui combinent la methode des elements finis et les methodes qui utilisent
les solutions fondamentales, par exemple la methode des elements de frontieres. Lavantage est
que la condition de radiation est deja automatiquement verifiee dans la solution elementaire.
Cependant, la solution nest pas toujours facile a trouver dans le cas general [18].
Ajouter une couche absorbante dans laquelle un amortissement assez important sera introduit
et de facon reguliere pour que les ondes puissent entrer dans cette couche sans etre reflechies et
que toutes les energies soient absorbees avant darriver a la frontiere [24][112].
Utiliser les procedures purement numeriques proposees par Wolf et Song [107]
etc.
Dans le cadre de lapplication aux problemes de charge mobile, la methode des elements aux
frontieres est construite de facon usuelle : Pan et al. [90], Aubry et al. [8] ... Recemment, Rassmusen
et al. [95] ont formule la fonction de Green pour le deplacement et la traction sur la surface dun
demi-espace dans le repere mobile. Les resultats pour le cas 2D sont compares a ceux dun calcul
en elements finis utilisant une condition absorbante qui est proposee par Krenk et al. [65].
Dans leur calcul, Krenk et al. ont introduit une matrice dimpedance mise a la frontiere en
tenant compte du terme convectif qui vient du fait que le probleme est passe dans le repere mobile.
30 1. Synthese bibliographique
Une matrice non-symetrique a ete obtenue qui est formulee a partir de deux directions separees pour
londe P et pour londe S. Une technique similaire appliquee au probleme dune source acoustique
mobile est aussi utilisee par Kirkegaard et al. dans [61].
Toujours dans le repere mobile, mais pour le probleme dune poutre dEuler-Bernouilli posee
sur la fondation de Winkler, Andersen et al. [7] ont propose une condition absorbante qui est basee
sur une relation explicite entre la force de reaction, le deplacement et la vitesse a la frontiere dans
le domaine du temps. Cette relation est obtenue par une linearisation de la fonction de transfert
force-deplacement dans le domaine des frequences. Cela donne le resulat exact si la charge est
harmonique mais pas avec un signal de charge quelconque.
Le ballast est un materiau granulaire. Dans la nature, les materiaux granulaires (mais aussi les
sols, poudres, roches, ...) sont des milieux biphasiques constitues dune phase solide dispersive et
dune phase fluide. La propriete essentielle qui distingue un milieu granulaire par rapport aux autres
milieux solides multiphasiques est la nature discrete de la phase solide. Comme la distribution
statique et dynamique des contraintes depend de letat des contacts entre les particules au moment
considere, le milieu granulaire est un milieu tres nonlineaire et dissipatif.
On va presenter dans la suite deux aspects pour modeliser le materiau granulaire : le modele
discret et le modele continu.
Ce modele considere le milieu ballast comme un ensemble delements discrets. Chaque element
peut etre deformable ou indeformable. Lhypothese de corps rigide est la plus utilisee puisque
dans une condition de contact normale, les grains dans le milieu granulaire sont tres sensibles aux
mouvements dans les directions de traction et de glissement.
La methode pour etudier numeriquement le probleme dynamique dans ce modele est la methode
des elements discrets (MED). En general, en utilisant lhypothese des particules rigides, MED part
des equations de mouvement de Newton ou chaque particule i dans lensemble de N particules a
6 degres de liberte (3 en translation et 3 en rotations) [35] :
N
X N
X
mi ui = f gi + f cij ; Ii i = mgi + mcij (1.23)
j=1 j=1
ou mi et Ii sont la masse et le moment inertie de la particule i, f cij et mcij sont les vecteurs
forces et moments de contacts appliques sur la particule, f gi et mgi sont les forces et les moments
autre que le contact (volumique) sur la particule i. Des algorithmes dintegrations numeriques
(methode de difference finie) sont utilises pour resoudre ce systeme. La difficulte est la detection
des contacts pour chaque particule a chaque pas de temps. On peut citer plusieurs etudes sur
le schema dintegration et lalgorithme de detection, par exemple : Azanza [9], Wassgren [106],
Oviedo [88][89] pour les elements circulaires 2D, Oviedo [88] pour les elements triangulaires 2D,
Cundall [25], Ghraboussi [42] pour les elements polyhedriques 3D ... Il a montre que, en prenant
un pas de temps de calcul suffisamment petit dans la simulation numerique, les ondes ne peuvent
pas se propager entre les grains voisins, et la vibration dune particule ne depend que des forces
de contact (voir [106] [113]).
La plupart des etudes realisees avec la methode des elements discrets sont construites pour
etudier la vibration dun ensemble de particules. Le probleme de propagation dondes pour evaluer
la reponse dynamique des voies ballastees est peu etudie. Dans la litterature, le probleme des ondes
1.6 Conclusion 31
dans linterieur des particules nest pas vraiment traite, mais en fait, la propagation des ondes dans
un paquet de grains est consideree comme la propagation des contraintes aux points de contacts.
Citons par exemple ; les travaux de Melin [79] (qui a montre que la vitesse des ondes depend aussi
de lamplitude de la force appliquee), de Sadd et al. [98] (qui ont montre en etudiant des particules
elliptiques que la propagation des ondes dans un ensemble de particules depend beaucoup des
arrangements spatiaux) , de Shulka et al. [102] (qui ont compare les resultats experimentaux et
numeriques sur un ensemble de particules circulaires) ...
Ce modele est realiste pour modeliser les materiaux granulaires mais par contre, demande un
volume de calcul tres important.
Dans le cas ou la taille des particules est suffisamment petite par rapport a la dimension du
domaine, le modele continu peut etre utilise pour etudier le probleme global. Les modeles continus
construisent les lois de comportement a partir de la microstructure discrete (equation 1.23). Granik
et Ferrari [43] ont propose un modele elastique lineaire, puis dans [78], Maddalena et Ferrari ont
propose un modele de viscoelaticite. Des modeles elastiques dordre plus eleve ont ete proposes par
Chang et al. [19][20], Muhlhaus [84] ... En outre, une discussion sur la definition des contraintes et
des deformations dans un domaine de materiaux granulaires a ete donnee par Bagi [10].
Le comportement dynamique du materiau granulaire a ete etudie en utilisant ces modeles conti-
nus. Les equations de dispersion dondes presentees dans [19] montrent que le materiau granulaire
est par nature un filtre dans lequel les ondes de petites longueurs et de hautes frequences ne peuvent
pas se progager. Le probleme dynamique dune couche de materiau granulaire posee sur un demi-
espace est considere dans les travaux de Suiker et al. [103] pour une charge fixe ou [104][105] pour
une charge mobile.
1.6 Conclusion
Bien que le probleme dynamique induit par le trafic a ete beaucoup considere dans la litterature,
tres peu detudes ont ete faites en tenant compte de la non-linearite de la structure. On va proposer
dans cette these, des modeles et des methodes de resolution qui permettent devaluer linfluence de
la non-linearite sur la reponse dynamique des structures due aux charges mobiles. Les valeurs des
vitesses critiques avec la presence de la non-linearite seront aussi determinees. Un modele continu
qui permet de tenir compte de leffet unilateral dans le ballast sera aussi introduit.
32 1. Synthese bibliographique
Chapitre 2
2.1 Introduction
2 u(x, t)
2 u(x, t) + ( + ) u(x, t) + g = (2.1)
t2
ou u(x, t) est le vecteur deplacement au point x a linstant t (u = {u 1 , u2 , u3 }t ). Les coefficients
, , sont respectivement la densite et les deux coefficients elastiques de Lame de la couche
consideree. Le terme g designe la force volumique qui sera omise dans la suite car seule la partie
dynamique de la reponse sera calculee.
34 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
repere fixe
repere mobile
x2 x1 f0 (t)
x3 vt v
X1
X2
X3
1 , E 1 , 1
2 , E 2 , 2
La charge mobile sapplique sur la surface et se deplace avec une vitesse constante v. Cest une
charge harmonique et elle peut avoir un des deux types suivants :
soit une force ponctuelle :
soit une force uniforme sur une surface ([a, a] [b, b]) :
Comme la charge se deplace avec une vitesse constante, la solution en regime permanent de
(2.1) devient stationnaire dans un repere mobile qui se deplace avec la meme vitesse. Pour passer
dans ce nouveau repere, on realise le changement de variable :
x X = x vte1 (2.5)
La solution harmonique stationnaire u(x, t) dans ce repere mobile peut etre representee par :
u U
= iU v (2.7)
t X1
2
u U 2U
= 2 U 2iv + v2 (2.8)
t2 X1 X12
En remplacant (2.6) et (2.8) dans (2.1), lequation elastodynamique dans le repere mobile devient
une equation differentielle statique :
U 2U
2 U + ( + ) U + 2 U + 2iv v 2 =0 (2.9)
X1 X12
U 2U
c2s 2 U + (c2p c2s ) U + 2 U + 2iv v2 =0 (2.10)
X1 X12
Comme le milieu considere est homogene dans deux directions horizontales (Oe 1 e2 ), on va faire
une transformation de Fourier des equations dans ces deux directions :
Z + Z +
ik2 X2
U (k1 , k2 , X3 ) = e dX2 U (X1 , X2 , X3 )eik1 X1 dX1 (2.12)
Z + Z +
1
U (X1 , X2 , X3 ) = eik2 X2 dk2 U (k1 , k2 , X3 )eik1 X1 dk1 (2.13)
(2)2
En remplacant (2.13) dans (2.10), on peut ecrire un systeme dequations differentielles sur U dans
le domaine des nombres donde par rapport a la seule variable X 3 :
2 U (k1 , k2 , X3 ) U (k1 , k2 , X3 )
A + iB C U (k1 , k2 , X3 ) = 0 (2.14)
X32 X3
ou :
c2s 0 0 0 0 k1 (c2p c2s )
A = 0 c2s 0 ; B = 0 0 k2 (c2p c2s ) (2.15)
0 0 c2p k1 (c2p c2s ) k2 (c2p c2s ) 0
k12 c2p + k22 c2s ( vk1 )2 k1 k2 (c2p c2s )
C = k1 k2 (c2p c2s ) k12 c2s + k22 c2p ( vk1 )2 0 (2.16)
0 0 2 2 2 2
k1 cs + k2 cs ( vk1 ) 2
36 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
Demonstration.
On evalue chaque partie de (2.10) dans le domaine des nombres donde :
2U 2U 2U TF 2 U
2 U = + + k 2
1 U k 2
2 U + (2.17)
X12 X22 X32 X32
2
U1 2U 2 2U 3 2 U 3
X 2 + + k1 U 1 k1 k2 U 2 + ik1 X
1 X 1 X 2 X 1 X 3
3
2U 1 2
U2 2
U3 TF
U 3
( U ) = 2
X X + X 2 + X X k1 k2 U 1 k2 U 2 + ik2 X
(2.18)
2 1 2 2 3 3
2U 1 2U 2 2U 3 U 1 U 2 2
U 3
+ + ik1 + ik2 +
X3 X1 X3 X2 X32 X3 X3 X32
U 2U T F 2
2 U + 2iv v2 U 2vk1 U + v 2 k12 U (2.19)
X1 X12
En faisant la somme de ces trois expressions et en mettant sous forme matricielle, on obtient (2.14).
On peut aussi exprimer la TF des trois composantes du tenseur des contraintes sur un plan hori-
zontal (t = .n3 ) :
t = {13 , 23 , 33 }t
( ! ! )t
U1 U2 U3
= ik1 U3 + , ik2 U3 + , ik1 U1 + ik2 U2 + ( + 2) (2.20)
X3 X3 X3
det(k32 A + k3 B + C) = 0 (2.22)
soit :
c2s k32 + k12 c2p + k22 c2s ( vk1 )2 k1 k2 (c2p c2s ) k1 k3 (c2p c2s )
det k1 k2 (c2p c2s ) c2s k32 + k12 c2s + k22 c2p ( vk1 )2 k2 k3 (c2p c2s ) =0
k1 k3 (c2p c2s ) k2 k3 (c2p c2s ) c2p k32 + k12 c2s + k22 c2s ( vk1 )2
(2.23)
Ce determinant se reduit a une expression assez simple apres quelques simplifications :
[c2p (k12 + k22 + k32 ) ( vk1 )2 ][c2s (k12 + k22 + k32 ) ( vk1 )2 ]2 = 0 (2.24)
ou :
v v
kp = ; ks = ; mp = ; ms = (2.27)
cp cs cp cs
2.3 Methode de resolution par la transformation de Fourier 37
q q
p = k12 (kp mp k1 )2 + k22 ; s = k12 (ks ms k1 )2 + k22 ; (2.28)
On a deux valeurs propres uniques pour londe de compression et deux valeurs propres doubles
pour les ondes de cisaillement. Sans perte de generalite, on suppose que : Re ( p ) > 0 ; Re (s ) > 0.
Les deux constantes mp et ms sont connues dans la litterature comme les nombres de Mach.
Les vecteurs propres correspondants a k 3 sont :
k1 0 s
U1,2 = k2 U3,4 = s U5,6 = 0
ip ik2 ik1
ip A+ 1e
p X3
+ ik2 A+2e
s X3
+ ik1 A+3e
s X3
(2.31)
et on peut aussi en deduire les expressions de .n3 :
13 = [2k1 p A
1e
p X3
k 1 k2 A
2e
s X3
(k12 + 2s )A
3e
s X3
(2.32)
+2k1 p A+ 1e
p X3
k 1 k2 A+
2e
s X3
(k12 + 2s )A+3e
s X3
]
p X3 2 2 s X3 s X3
23 = [2k2 p A1 e (k2 + s )A2 e k 1 k2 A3 e (2.33)
+ p X3 + s X3
+2k2 p A1 e 2 2
(k2 + s )A2 e k 1 k2 A+3e
s X3
]
2 2 2 p X3 s X3
33 = i[(k1 + k2 + s )A1 e 2k2 s A2 e 2k1 s A 3e
s X3
(2.34)
(k12 + k22 + 2s )A+1e
p X3
+ 2k2 s A+ 2e
s X3
+ 2k1 s A+3e
s X3
]
Demonstration.
Deux expressions (2.32) et (2.33) sont evidentes a partir de (2.20) en derivant la solution generale de
U . La solution 33 vient de :
U3
33 = (ik1 U1 + ik2 U2 ) + ( + 2)
X3
= i [(k12 + k22 ) 2p ( + 2)]A
1 e
p X3
+ i [(k12 + k22 ) 2p ( + 2)]A+
1e
p X 3
+i[2k2 s A
2e
s X3
+ 2k2 s A+
2e
s X 3
2k1 s A
3e
s X3
+ 2k1 s A+
3e
s X 3
] (2.35)
On va simplifier encore le terme (k12 + k22 ) 2p ( + 2) de la facon suivante :
2 2 2 2 2 2 2 2 2 2 ( vk1 )2
(k1 + k2 ) p ( + 2) = (cp 2cs )(k1 + k2 ) cp k1 + k2
c2p
= 2c2s (k12 + k22 ) + ( vk1 )2
= c2s (k12 + k22 ) c2s (k12 + k22 ) ( vk1 )2
= [k12 + k22 + 2s ] (2.36)
On en deduit (2.34)
Pour la commodite du calcul, on va representer les expressions de U et t sous la forme matricielle
suivante :
U
= CeA
t
C11 C12 e 0 A
= (2.37)
C21 C22 0 e+ A+
38 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
ou :
A = {A t
1 , A2 , A3 } ; A+ = {A+ + + t
1 , A2 , A3 } ; (2.38)
ep X3 0 0 e p X3 0 0
e = 0 es X3 0 ; e+ = 0 e s X3 0 (2.39)
0 0 e s X3
0 0 e s X3
k1 0 s k1 0 s
C11 = k2 s 0 C12 = k2 s 0 (2.40)
ip ik2 ik1 ip ik2 ik1
2k1 p k1 k2 (k12 + 2s )
C21 = 2k2 p (k22 + 2s ) k1 k2 (2.41)
2 2 2
i(k1 + k2 + s ) 2ik2 s 2ik1 s
2k1 p k1 k2 (k12 + 2s )
C22 = 2k2 p (k22 + 2s ) k1 k2 (2.42)
2 2 2
i(k1 + k2 + s ) 2ik2 s 2ik1 s
ou f est la TF du vecteur force applique sur la surface qui peut etre determinee en fonction de
sa forme :
2.3 Methode de resolution par la transformation de Fourier 39
Condition aux limites a linfini. La derniere couche n est un demi-espace et donc quand
X3 + : U 0. Cette condition impose que :
A+
n =0 (2.48)
ou A est toujours constant, le probleme ne change donc pas si lon choisit localement pour chaque
couche i une origine des coordonnees a linterface en haut [i 1, i]. Lequation de la condition de
continuite entre deux couches [i, i + 1] (2.44) ecrite avec ce choix devient :
i i+1
C11 Ci12 ei 0 Ai C11 Ci+1 12 Ai+1
= (2.50)
Ci21 Ci22 0 e+ i A +
i C i+1
21 C i+1
22 A+ i+1
ou : e
i = ei (X3 = di ) et di est lepaisseur de la couche i. On note aussi que comme cette condition
est ecrite a lorigine des coordonnees de la couche i + 1, la matrice exponentielle est eliminee a
+
droite puisque : e i+1 (X3 = 0) = ei+1 (X3 = 0) = 1.
Il est plus interessant decrire la condition de continuite (2.50) comme une relation entre les
ondes entrantes et les ondes sortantes de linterface [i, i + 1] comme suit :
A+i i,i+1 Ai
=Q (2.51)
A
i+1 A+
i+1
A
i1 A+
i1
couche i-1
A
i A+
i
couche i
couche i
A+ A
i+1 couche i+1
i+1
A+
i+1 A
i+1 couche i+1
ou :
i+1 1 i
i+1 1 i 1
Li,i+1
11 = (Ci+1 1 i
11 ) C12 (C21 ) C22 (C11 ) C11 (Ci+1 1 i
21 ) C21 (2.53)
i+1 1 i+1
i+1 1 i 1
Li,i+1
12 = (Ci+1 1 i
11 ) C12 (C21 ) C22 (C11 ) C12 (Ci+1 1 i+1
21 ) C22 (2.54)
i 1 i
i 1 i+1 1
Li,i+1
21 = (Ci12 )1 Ci+1
11 (C22 ) C21 (C12 ) C11 (Ci22 )1 Ci21 (2.55)
i 1 i+1
i 1 i+1 1
Li,i+1
22 = (Ci12 )1 Ci+1
11 (C22 ) C21 (C12 ) C12 (Ci22 )1 Ci+1 22 (2.56)
Avec cette facon decrire la condition de continuite, on nest plus gene par les exponentielles
positives. Si lepaisseur de la couche i est importante, le terme e
i devient tres petit et cela conduit
a une tres petite onde reflechie dans cette couche. On peut aussi verifier que si deux couches sont
identiques, Ci = Ci+1 implique que L11 = 0, L12 = 1, L21 = 1, L22 = 0, i.e. A+ +
i = ei Ai+1 et
+
A i = ei Ai+1 .
Par consequent, en connaissant les relations (2.51) pour deux interfaces [i 1, i] et [i, i + 1], on
peut determiner la relation entre les ondes entrantes et les ondes sortantes dans la couche i :
A+ A
i1 = Qi1,i+1 i1 (2.57)
A
i+1 A+
i+1
Qi1,i+1
11 = Qi1,i
11 + Qi1,i
12 Ria Qi,i+1
11 Qi1,i
21 (2.59)
Qi1,i+1
12 = Qi1,i
12 Ria Qi,i+1
12 (2.60)
Qi1,i+1
21 = Qi,i+1
21 R i
Q
b 21
i1,i
(2.61)
Qi1,i+1
11 = Q21 Rb Qi1,i
i,i+1 i
22 Qi,i+1
12 + Qi,i+1
22 (2.62)
1 1
ou Ria = 1 Qi,i+1
11 Q22
i1,i
et Rib = 1 Qi1,i
22 Q11
i,i+1
.
2.3 Methode de resolution par la transformation de Fourier 41
Demonstration.
Pour simplifier les expressions, on va noter Qi1,i et Qi,i+1 respectivement par Q et Q+ . En utilisant
lequation (2.51) :
+
Ai1 Ai1 A+i + A i
=Q et =Q (2.63)
Ai A+i Ai+1 A+i+1
Lorsque A + + +
i et Ai sont connus, on peut exprimer Ai1 et Ai+1 en fonction de Ai1 et Ai+1 en remplacant
+
Ai , Ai dans (2.63) :
A+
i1 = Q +
11 Ai1 + Q12 Ai
+
1
= Q +
11 Ai1 + Q12 1 Q11 Q22 Q11 Q21 Ai1 + Q+ +
12 Ai+1
h i
1 1
= Q +
11 + Q12 1 Q11 Q22 Q+ +
11 Q21 Ai1 + Q12 1 Q11 Q22 Q+ +
12 Ai+1 (2.65)
A
i+1 = Q+ + +
21 Ai + Q22 Ai+1
+ 1
= Q+
21 1 Q22 Q11 Q21 Ai1 + Q + + + +
22 Q12 Ai+1 + Q22 Ai+1
h i
+ 1 + 1
= Q+
21 1 Q22 Q11 Q +
21 Ai1 + Q21 1 Q22 Q11 Q + + +
22 Q12 + Q22 Ai+1 (2.66)
1 1
En notant Ria = 1 Q+
11 Q22 et Rib = 1 Q +
22 Q11 , on deduit (2.57).
En utilisant (2.51) et (2.57), on est capable de calculer la matrice de transfert Q pour un
nombre de couches quelconque de facon iterative. On peut presenter la relation entre les ondes
dans la premiere couche 1 et la derniere n par :
+
A1 1,n A1
=Q (2.67)
A n A+
n
Lorsque A1 est connu, les vecteurs Ai de toutes les couches peuvent etre calcules grace aux deux
relations (2.51) et (2.57). Supposons que A i soit deja connu, Ai+1 peut etre determine de la facon
suivante :
42 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
i,i+1 i,i+1 +
(2.51) A
i+1 = Q21 Ai + Q22 Ai+1
i+1,n
(2.57) A+
i+1 = Q11 Ai+1
donc :
1
i,i+1 i+1,n
A
i+1 = 1 Q 22 Q 11 Qi,i+1
21 Ai
1 (2.71)
i+1,n i,i+1 i+1,n
A+
i+1 = Q 11 1 Q 22 Q 11 Qi,i+1
21 Ai
Remarque 2.3.1
1. Il faut faire attention de ne pas calculer A +
i+1 et Ai+1 en utilisant directement lequation (2.51) parce
quelle va faire apparatre des termes exponentiels positifs. La formule (2.71) est plus compliquee mais
ne concerne que des exponentielles negatives.
2. Si le milieu considere est un demi-espace homogene (une seule couche), on a tout de suite A +
1 = 0 et
donc A1 est evalue simplement par :
1
A 1
1 = C21 f (2.72)
v = i( vk1 )U (2.75)
2
a = ( 2vk1 + v 2 k12 )U (2.76)
En connaissant la solution dans le domaine des nombres dondes, il faudra calculer la trans-
formation de Fourier inverse (2.13) pour obtenir la solution en espace. Afin de pouvoir evaluer
precisement cette integration de la transformation inverse, il est interessant de voir dabord quel
est le comportement de la solution dans le domaine des nombres donde et surtout de chercher sil
existe des points singuliers. On va considerer ici le cas plus simple ou le domaine est un demi-espace
homogene.
on cherche sil existe des points singuliers de cette fonction. Pour cela, il suffit de trouver les racines
du determinant de C21 qui sexprime par :
on rappelle que :
q q
p = k12 (kp mp k1 )2 + k22 ; s = k12 (ks ms k1 )2 + k22
v
kp,s = ; mp,s =
cp,s cp,s
On peut trouver tout de suite que k1 = k2 = 0 est une racine de (2.78) si et seulement si = 0.
Pour les autres racines, en supposant que k 12 + k22 6= 0, on peut diviser (2.78) par (k12 + k22 )2 pour
obtenir :
2 21 21
(ks ms k1 )2 (ks ms k1 )2 (kp mp k1 )2
2 4 1 1 =0 (2.79)
k12 + k22 k12 + k22 k12 + k22
soit :
!2 2 !2 12 !2 12
2 1 vk1
p 4 1 1 vk1
p 1 1 vk1
p =0 (2.80)
c2s 2
k1 + k22 c2s 2
k1 + k22 c2p k12 + k22
Cette equation est bien connue comme lequation de Rayleigh et a toujours une racine differente
de zero. La racine non-nulle cR est appelee la vitesse donde de surface (ou la vitesse donde de
Rayleigh). Alors, les points singuliers vont apparatre au lieu des points qui verifient :
| vk1 |
p = cR (2.81)
k12 + k22
Rappelons que k1 et k2 sont reels. Donc, la condition dexistence de k 2 qui verifie (2.81) implique :
v cR (2.83)
k1
Par ailleurs, les points qui verifient s = 0 sont aussi les racines de D(k1 , k2 ) (2.78). Ces points
doivent etre les racines dune equation qui a la meme forme que lequation 2.81 :
| vk1 |
p = cs (2.84)
k12 + k22
La condition dexistence de k2 qui verifie cette equation est : v cs .
k1
En bref, U(k1 , k2 ) nest pas du tout une fonction reguliere. Elle a des points singuliers appa-
raissant pour plusieurs valeurs de k 1 , k2 :
sur la ligne k1 = /v, k2 .
ensuite, il y a des possibilites davoir dautres points singuliers qui dependent des valeurs de
et de v :
si 6= 0, v < cR < cs : il y a des points singuliers sur deux couples de lignes qui sont limitees
dans k1 [/(cR v), /(cR + v)] et k1 [/(cs v), /(cs + v)] .
44 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
si 6= 0, cs > v > cR : il y a des points singuliers sur deux couples de lignes qui sont
respectivement definies telles que k 1 [, /(cR + v)] [/(v cR ), +] et k1 [/(cs
v), /(cs + v)].
si 6= 0, v > cs > cR : il y a des points singuliers sur deux couples de lignes qui sont
respectivement definies telles que k 1 [, /(cR +v)][/(vcR ), +] et k1 [, /(cs +
v)] [/(v cs ), +]
Une representation des positions des points singuliers dans deux cas subsoniques et supersoniques
est presentee dans la figure 2.3.
p+
s p+
s
v
p p p+ p p+ p
k2
k2
R s R v s R R
k1 k1
p
R = ; p+
R =
cR v cR + v
p
s = ; p+
s =
cs v cs + v
p+
s
v
p+
k2
R p
s
p
R
k1
Remarque 2.3.2
On peut en deduire facilement que dans le cas de la charge constante, i.e = 0, les points singuliers se
trouvent :
sur la ligne k1 = 0.
si v < cR < cs : il ny a pas dautres points singuliers.
p
(v/cR )2 1.
si cs > v > cR : il y a des points singuliers sur les deux lignes droites k2 = k1
p
si v > csp> cR : il y a des points singuliers sur les quatre lignes droites k2 = k1 (v/cR )2 1 et
k2 = k1 (v/cs )2 1.
2.3 Methode de resolution par la transformation de Fourier 45
integration de Gauss
k2
k1 k2
k1
Figure 2.4.
Aux points qui se trouvent sur la ligne k 1 = /v, la fonction h(/v, k2 ) est calculee de facon
approximative par :
1
h(/v, k2 ) = [h(/v , k2 ) + h(/v + , k2 )] (2.85)
2
ou designe la viscosite. Ce coefficient de viscosite deplace les poles de laxe reel et la fonction
dans lintegration de la transformation est definie numeriquement de facon normale.
Remarque 2.3.3
1. Dans le cas subsonique, les fonctions de U et sont assez regulieres, lalgorithme de transformation de
Fourier rapide (FFT) dans deux dimensions (k1 , k2 ) peut etre utilise avec un petit nombre de points.
Lutilisation de la FFT dans ce cas est preferee comme elle est beaucoup plus rapide que celle qui utilise
lintegration numerique.
46 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
2. On peut reduire le temps de calcul de lintegration numerique de Filon dans la direction k 2 en remarquant
que les composantes u11 , u13 , u22 , u31 , 13
1 2
, 23 3
, 33 3
, 13 sont symetriques en fonction du nombre donde k2
i
(on note (.) la reponse due a la force qui sapplique suivante la direction i).
3. Le module delasticite complexe secrit habituellement E = E (1 + i) pour les equations vibratoires
classiques dans le domaine des frequences. Cependant, le module dans ce cas (2.86) depend aussi du signe
de vk1 . En fait, ce terme vient de la loi de comportement visco-elastique ecrite dans le domaine des
nombres dondes apres un changement de variable. Considerons une relation de contrainte-deformation
visco-elastique uniaxiale du type de Kelvin-Voigt suivante :
= E + (2.87)
Cette relation dans le repere mobile en remplacant par v,x secrit : = E + (i v,x ) et sa
transformation dans le domaine des nombres donde donne :
= E
+ i( vk1 )
(2.88)
on voit bien que le module delasticite est complexe et depend strictement du signe de ( vk 1 ).
Le premier test est pour valider la methode de calcul dans le cas plus simple ou on considere une
force statique fixee sappliquant sur la surface libre dun demi-espace homogene. Ceci est appele
0.0002 0.0018
numerique numerique
analytique analytique
0.00015 0.0016
0.0014
1e-04
0.0012
Deplacement u1 (m)
Deplacement u3 (m)
5e-05
0.001
0
0.0008
-5e-05
0.0006
-0.0001
0.0004
-0.00015 0.0002
-0.0002 0
-4 -2 0 2 4 -4 -2 0 2 4
X1 (m) X1 (m)
2000 0
numerique numerique
analytique analytique
1500
-1000
1000
-2000
Contraintes 13 (N/m2)
Contraintes 33 (N/m2)
500
0 -3000
-500
-4000
-1000
-5000
-1500
-2000 -6000
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5
X1 (m) X1 (m)
le probleme de Boussinesq et ses solutions sont donnees explicitement sous forme analytique (voir
par exemple dans [56], chapitre 3).
Les caracteristiques physiques du milieu sont la masse volumique de 1000 kg.m 3 , le module
dYoung E de 50 MPa et le coefficient de Poison de 0.3. La force statique a une amplitude de
1000 N et applique verticalement sur la surface libre du milieu solide. A cause de la singularite qui
apparat quand la vitesse tend vers zero, la force statique est simulee par une charge qui se deplace
avec une tres petite vitesse de 0.5 m.s 1 .
Les nombres dechantillons dans les deux directions k 1 , k2 sont pris identiques N1 = N2 = 256
points et les pas dechantillonage sont k 1 = k2 = 0.2. Pour lintegration numerique au voisinage
de la force, on prend 7 points de Gauss.
On presente dans la figure 2.5 la comparaison entre les resultats numeriques et analytiques.
Les solutions comparees sont des deplacements (horizontaux et verticaux) et des contraintes (de
cisaillement et normales) a la profondeur de h = 0.3 m. Les resultats obtenus sont parfaitement
identiques dans les deux cas.
On profite aussi de cet exemple pour voir leffet de lintegration par les points de Gauss au
voisinage du point singulier k1 = k2 = 0. La difference entre deux calculs avec ou sans lintegration
de Gauss est presentee dans la figure 2.6.
0.0018
avec pts de Gauss
sans pts de Gauss
analytique
0.0016
0.0014
0.0012
Deplacement u3 (m)
0.001
0.0008
0.0006
0.0004
0.0002
0
-4 -2 0 2 4
X1 (m)
On va considerer dans cet exemple le cas dun demi-espace homogene. Les caracteristiques
physiques de ce demi-espace sont presentees dans le tableau 2.1. Le probleme est calcule dabord
avec une charge constante, puis avec une charge harmonique.
A. Charge constante
On applique une charge mobile dont lamplitude est de 170 kN a differentes vitesses. En com-
parant aux vitesses des ondes de Rayleigh, de cisaillement et de compression (tableau 2.1), on va
utiliser quatre valeurs de la vitesse de charge : 100 m.s 1 , 150 m.s1 , 200 m.s1 et 300 m.s1 afin
dobtenir les cas subsonique, transonique et supersonique (ces quatre vitesses correspondent aux
cas v < cR , cR < v < cs , cs < v < cp et v > cp ).
8 25
v=100 m/s v=100 m/s
6 v=150 m/s v=150 m/s
v=200 m/s v=200 m/s
v=300 m/s 20 v=300 m/s
4
2 15
Deplacement u1 (mm)
Deplacement u3 (mm)
0
10
-2
5
-4
-6 0
-8
-5
-10
-12 -10
-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4
X1 (m) X1 (m)
50 40
v=100 m/s v=100 m/s
40 v=150 m/s 20 v=150 m/s
v=200 m/s v=200 m/s
v=300 m/s v=300 m/s
30 0
20 -20
Contraintes 13 (MPa)
Contraintes 33 (MPa)
10 -40
0 -60
-10 -80
-20 -100
-30 -120
-40 -140
-50 -160
-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2
X1 (m) X1 (m)
La figure 2.7 presente les resultats des deplacements horizontaux et verticaux, aussi des contraintes
normales et de cisaillement dus aux differentes vitesses. Le profondeur de calcul est 0.1 m. Les pa-
rametres de calcul comme le nombre des points pour la FFT, la frequence cut-off des nombres
donde ... sont choisis en fonction des valeurs des vitesses. Le coefficient de viscosite utilise dans
les cas transonique ou supersonique est de 2.5%.
On constate que les points maximums sont deplaces vers la gauche de la position de la charge
quand la vitesse augmente. Les valeurs de deplacement et de contrainte sont plus elevees mais
gardent toujours les memes formes par rapport a celles du cas statique. En revanche, des que la
vitesse de la charge depasse la vitesse donde de surface, les formes sont vraiment changees. La
reponse devant la force est de plus en plus diminuee et dans le cas supersonique (v = 300 m.s 1),
il ny a plus rien dans cette partie.
La reponse est plus importante dans les cas transoniques, surtout quand la vitesse de la charge
2.4 Validation et exemples numeriques 49
est au voisinage de la vitesse donde de Rayleigh. Quand v = 200 m.s 1 , cs < v < cp , la reponse
devant la force nest due quaux ondes de compression qui se propagent radialement. Cest pourquoi
les resultats de ce cas montrent quil y a presque uniquement le deplacement horizontal dans la
zone X1 > 0 .
v = 100 m/s
u3 (mm)
14
12
10
8
6
4
2
0
6
4
2
8 0
6
4 2 Y (m)
2
0 4
2
X (m) 4
6 6
v = 300 m/s
u3 (mm)
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0.2
0.4
8
6
4
2
20 0
15
10 2 Y (m)
5
0 4
5
10 6
X (m) 15
20 8
La figure 2.9 presente les resultats dacceleration verticale. On peut constater que les valeurs
maximums des accelerations deviennent tres grandes (impossible dans la realite) dans les cas su-
personiques aux endroits ou il y a des chocs. Ce probleme vient de la singularite due a la force
ponctuelle. On va voir apres que les resultats peuvent etre nettement ameliores en utilisant une
charge de forme reguliere (de type Gaussien).
50 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
10000 1.5e+06
v = 100m/s v = 150m/s
5000
1e+06
0
-5000
Acceleration a3 (m.s-2)
Acceleration a3 (m.s-2)
500000
-10000
-15000 0
-20000
-500000
-25000
-30000
-1e+06
-35000
-40000 -1.5e+06
-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4
X1 (m) X1 (m)
1.5e+07 1e+07
v = 200m/s v = 300m/s
1e+07
5e+06
5e+06
Acceleration a3 (m.s-2)
Acceleration a3 (m.s-2)
0 0
-5e+06
-5e+06
-1e+07
-1.5e+07 -1e+07
-2e+07
-1.5e+07
-2.5e+07
-3e+07 -2e+07
-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4
X1 (m) X1 (m)
B. Charge harmonique
Le probleme est considere cette fois avec une charge mobile ponctuelle harmonique dont la
frequence est f0 (f0 = ). Les calculs sont realises avec differentes vitesses et frequences de
2
la charge pour voir ses influences sur le domaine. On presente ici seulement les solutions de
deplacement vertical.
Effet des vitesses. Une force harmonique de frequence 20Hz est utilisee dans le calcul pour les
differentes vitesses de v = 50, 100, 150, 200, 300m.s 1 , les deplacements sont calcules a la profon-
deur 0.25m.
Les solutions quand v = 50 et 100m.s1 presentees dans la figure 2.10 montrent bien leffet
Doppler ou la vitesse de la charge change la frequence de la propagation donde dans les deux
parties : devant et derriere la force. Les longueurs donde changent strictement en fonction de la
vitesse de charge.
Effet des frequences. On presente les resultats a la profondeur 0.05m pour plusieurs valeurs de
frequences (f0 = 20, 50, 100, 200 Hz) dans deux cas : subsonique (v = 100ms 1 ) et supersonique
(v = 300m.s1 ).
Les figures 2.11a,b montrent les resultats dans les echelles normales. Pour comparer les valeurs
aux points maximums, on trace dans les figures 2.11c,d les deplacements en fonction dun parametre
X1 f0
sans dimension defini par qui permet de placer les pics aux memes positions.
v
2.4 Validation et exemples numeriques 51
5 12
v = 50m/s v = 150m/s
v = 100m/s
5 12
10
4 4 10
3 8
6
2 8
4
3
Deplacement u3 (mm)
Deplacement u3 (mm)
1 2
0 0
6
-1 -2
2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2
1
2
0
0
-1 -2
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
X1 (m) X1 (m)
14 6
v = 200m/s v = 300m/s
12 14 6
12 5
5
10 4
10 8
3
6 4
4 2
Deplacement u3 (mm)
Deplacement u3 (mm)
8 2 1
0 3 0
-2 -1
6 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2
2
4
1
2
0 0
-2 -1
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
X1 (m) X1 (m)
Figure 2.10. Comparaison des reponses dun demi-espace homogene dues aux differentes vitesses.
6
8
4
u3 (mm)
u3 (mm)
6
2
4
0
2
-2
-4 0
-6 -2
-1.5 -1 -0.5 0 0.5 -1.5 -1 -0.5 0 0.5
X1 (m) X1 (m)
-1 -1
(c) v = 100m.s (d) v = 300m.s
12 14
12 f = 20Hz 14 f = 20Hz
10 f = 50Hz 12 f = 50Hz
10 f = 100Hz 12 f = 100Hz
8 f = 200Hz 10 f = 200Hz
6 8
8
4 10
6
2
6 0 4
8 2
-2
4 -4 0
u3 (mm)
u3 (mm)
-6 6 -2
-0.2 -0.1 0 0.1 0.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2
2
4
0
2
-2
-4 0
-6 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 -2.5 -2 -1.5 -1 -0.5 0 0.5 1
X1f0/v X1f0/v
Cet exemple considere le cas ou le milieu se compose de 2 couches de deux materiaux differents :
la premiere est une couche depaisseur 30cm, la deuxieme est un demi-espace infini. Les deux
materiaux utilises ici nont que les modules dYoung de differents. On note E 1 le module dYoung
de la premiere couche et E2 celui du demi-espace. On va considerer deux cas differents pour la
relativite entre les rigidites de deux couches :
(i) Une couche souple sur un demi-espace rigide : E 1 = 100 MPa, E2 = 200 MPa.
(ii) Une couche rigide posee sur un demi-espace souple : E 1 = 200 MPa, E2 = 100 MPa
Les caracteristiques qui correspondent a ces deux materiaux sont presentees dans le tableau 2.2.
Deplacement u3 (mm)
-5
8 10 -10
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
6 5
4 0
2 -5
0 -10
-1 -0.5 0 0.5 1 -3 -2 -1 0 1 2
X1 (m) X1 (m)
Deplacement u3 (mm)
3
8 0 2
5 1
-2
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0
6 4 -1
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1
3
4
2
2
1
0
0
-2 -1
-3 -2 -1 0 1 2 -3 -2 -1 0 1 2
X1 (m) X1 (m)
La figure 2.12 presente les solutions dans les deux cas en les comparant avec celles des deux
cas de milieu homogene : E = 100 MPa et E = 200 MPa. Les calculs sont realises pour 4 vitesses
differentes : v = 100, 200, 300, 400 m.s 1 . Le parametre evalue est le deplacement vertical dans la
premiere couche a la profondeur 0.1 m.
2.5 Applications 53
On appelle les solutions des cas (i) et (ii) respectivement u 1 et u2 et les solutions dans les cas
de milieux homogenes (E = 100, 200 MPa) sont u 1 et u2 . En regardant les resultats obtenus, on
peut faire les remarques suivantes :
- Quand v = 100 m.s1 , cette vitesse est subsonique pour les deux milieux. On constate que
u1 < u1 et u2 > u2 au voisinage de la force mais u1 = u1 et u2 = u2 au loin.
- Quand v = 200 m.s1 :
Pour le cas (i) : c1s < v < c2s , u1 = u1 au voisinage de la force.
Pour le cas (ii) : c1s > v > c2s , u2 et u2 ont les memes formes mais avec differentes valeurs.
- Quand v = 300 m.s1 ou v = 400 m.s1 : dans les deux cas, on a u1 = u1 et u2 = u2 au voisinage
de la force. Les differences deviennent plus importantes au loin.
En bref, on peut dire que la reponse au voisinage de la force dans la premiere couche est un
peu influencee par le demi-espace en dessous dans le cas subsonique. Au contraire, dans le cas
supersonique, elle est presque identique par rapport a celle du cas homogene. Physiquement, cest
parce que dans le cas supersonique, les ondes reflechies dues a la deuxieme couche nont pas
suffisamment de temps pour influencer les points au voisinage de la force.
2.5 Applications
On va traiter dans cette partie un probleme avec le chargement dun boggie qui se compose de
deux essieux comme dans la figure 2.13. Linteret de ce calcul est de trouver les cas critiques ou
les accelerations sur la surface de la couche de ballast depassent la gravite ce qui peut faire sauter
les cailloux de la voie.
Q = 17 t Q = 17 t
2b
Figure 2.13.
Pour calculer la force, Alaoui et Naciri [4] ont propose une formule analytique de la force
verticale exercee par le rail sur un blochet (i.e. la force totale dans le rectangle pointille dans la
figure 2.13) lors du passage dun boggie :
QY h ( V ta )2 V taL 2
i
F (t, Q, V ) = X d + X( d ) (2.89)
2
ou : Q est la charge par essieu (N ) qui vaut 17 tonnes, V est la vitesse du train (m.s 1 ), t est le
temps (s), d est la distance entre les traverses (d = 0.6m), a est une distance critique (a = 3m qui
est equivalent a 5d), X, Y sont des variables sans dimension comprises entre 0 et 1 qui dependent
du module dYoung Es du sol. Les valeurs de X et Y dependantes du module dYoung du sol sont
presentees dans le tableau 2.3.
54 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
Supposons que la charge est repartie uniformement sur les traverses dont la largeur est 2b =
QY 1
2.415m, la charge en espace a chaque instant t peut etre representee par :
V tx a 2 V tx1 aL 2
f (x1 , x2 , t) = X d
+X d
H(b |x2 |) (2.90)
2b
cette charge dans le repere mobile secrit :
! !
" "
2 2
X1 + L
2 X1 L
2
QY d d
f (X1 , X2 ) = X +X H(b |X2 |) (2.91)
2b
La figure (2.14) montre les signaux de force en temps (2.89) et en espace (2.91) qui sont sous la
forme dune lettre M .
40000 35000
35000
30000
30000
25000
25000
20000
force (N/m2)
force (N/m)
20000
15000
15000
10000
10000
5000
5000
0 0
0 0.05 0.1 0.15 0.2 -4 -2 0 2 4
temps (s) X1 (m)
Le couplage entre le sol et le ballast est considere comme une structure ou il y a une couche
depaisseur de 0.3m (pour le ballast) posee sur un demi-espace homogene (pour le sol). Cela conduit
a un probleme avec une distribution de charge mobile (Eq. 2.91) qui se deplace sur la surface dun
milieu a 2 couches. Il faut alors exprimer la formule analytique de la transformation de Fourier
(2.12) du signal de la charge (2.91) :
r 2 2
d k1 k1 L sin(k2 b)
f (k1 , k2 ) = 2QY d exp cos (2.92)
ln X 4 ln X 2 k2 b
On va presenter dans la suite seulement les solutions dacceleration verticale. Rappelons que la
direction positive de laxe vertical X 3 est vers en bas, la condition pour que les cailloux ne sautent
pas est amin > g ou g = 9.81 est lacceleration de gravite.
Les figures 2.15a et 2.15b montrent que les accelerations augmentent rapidement en fonction
de la vitesse et aussi de la frequence du train. Avec la presence de la frequence, les accelerations
ne sont plus symetriques par rapport a la position de la force et linfluence de la frequence est plus
importante sur la partie droite de la charge.
On remarque aussi que dans ce cas ou le signal de force est une fonction reguliere, il ny a plus
de probleme de points singuliers dus a la charge ponctuelle dans la solution en acceleration (comme
on a vu dans la figure 2.9).
5
2
0
Acceleration (m.s-2)
Acceleration (m.s-2)
0
-2
-4
-5
-6
-8
-10
-1
v = 75 ms -10
v = 100 ms-1 f0 = 5 Hz
v = 125 ms-1 f0 = 10 Hz
-15 -12
-15 -10 -5 0 5 10 15 -30 -20 -10 0 10 20 30
X1 (m) X1 (m)
On realise plusieurs calculs avec differentes caracteristiques physiques pour le ballast et le sol.
La charge est constante. On choisit deux modules dYoung pour le ballast (E b = 100, 200 MPa) et
trois modules dYoung pour le sol (E s = 60, 100, 150 MPa).
-2 -2
-4
-4
Acceleration min. (m.s-2)
-6
-6
-8
-8
-10
-10
-12
-12
-14
-14
-16
Pour chaque cas de (Eb ,Es ), on peut determiner les valeurs dacceleration minimum en fonction
de la vitesse de la charge. Ces fonctions sont tracees dans la figure 2.16 qui montrent des croissances
de type exponentiel des accelerations minimums quand la vitesse devient plus grande.
56 2. Charge harmonique mobile sur un demi-espace visco- elastique stratifie
Les resultats obtenus montrent que, avec le meme module dYoung du ballast E b , lacceleration
minimum est plus importante si le sol est plus souple. Ceci verifie le fait que le phenomene de guide
donde est plus important si la rigidite relative entre le ballast et le sol est plus grande (il y a plus
dondes flechies dans ce cas). On voit bien aussi que les accelerations minimums sont plus faibles
dans le cas ou la raideur du ballast est petite avec le meme sol.
2.6 Conclusion
Une methode semi-analytique pour resoudre le probleme lineaire de la propagation des ondes
dans un milieu solide multicouche dues aux forces appliquees sur la surface en regime permanent
est presentee dans cette partie. Les charges constantes ou harmoniques sont possibles. Un code de
calcul a ete aussi programme en langage C++.
Les exemples numeriques illustrent les reponses de differents milieux qui dependent fortement
de la vitesse et de la frequence de la charge, et aussi de lelasticite du milieu. Ils ont indique aussi
que la charge ponctuelle nest pas bonne pour decrire la charge reelle du train a cause de sa forme
singuliere. Il faudra utiliser un signal regulier pour eviter ce probleme. Un signal de force en M a
ete utilise pour tenir compte de la repartion de force du rail sur les blochets.
Pour un calcul plus detaille, il faut poser un probleme de couplage du milieu multicouche avec
le rail (qui peut etre modelise par une poutre dEuler-Bernouilli ou de Timoshenko). En plus,
la periodicite causee par les blochets peut avoir un drole important et doit etre ajoutee dans le
modele. La methode semi-analytique presentee dans cette partie est toujours valable.
Cependant, cette methode est limitee aux problemes lineaires. Cest pourquoi, on napprofon-
dira pas letude de ce cas et on ne cherchera pas a inclure le rail et les blochets. Dans le cas
non-lineaire, il faudra introduire des methodes numeriques, par exemple, la methode des elements
finis. On va discuter ce probleme dans les chapitres suivants.
Deuxieme partie
Modeles unidimensionnels
non-lineaires
Chapitre 3
3.1 Introduction
Dans les etudes sur le probleme dynamique des voies ballastees soumises aux passages des trains,
la structure de couplage entre le ballast et le sol peut etre consideree par un modele simplifie ou
la voie ballastee est modelisee par une poutre nonlineaire et la fondation (le sol) est modelisee
par un systeme de ressorts-amortisseurs (la fondation de Winkler). La figure 3.1 represente une
illustration de ce modele simplifie avec une force mobile de vitesse constante v.
y
f (x vt)
x , S, E +
Figure 3.1. Modele dune poutre soumise a une force mobile sur une fondation de Winkler
y
f0 (x vt)
x , S, E +
(k, )
On considere une barre soumise a une force axiale qui se deplace avec une vitesse v constante.
Cette barre est tenue par un systeme uniforme de ressorts - amortisseurs (figure 3.2) dont la rigidite
est k > 0 et lamortissement est > 0. Lequation dequilibre dynamique en petite deformation a
un point x en fonction du temps t secrit comme suit :
= Ec si 0 (3.3)
= Et si >0 (3.4)
On remarque que la condition initiale nest plus necessaire pour le probleme stationnaire. Pour
simplifier la notation, les sont supprimes dans les expressions suivantes. Le probleme se ramene
donc a resoudre lequation differentielle stationnaire du 2 e ordre avec les conditions aux limites
suivantes :
u(x) u(x)
c,t v + ku(x) = f0 (x)
x x x (3.6)
u() = u(+) = 0
ou c,t = S(Ec,t v 2 ), v = v.
On note que la variable peut aussi etre representee par c,t = S(c2c,t v 2 ), ou cc,t sont les
celerites des ondes de compression ou de traction :
s
Ec,t
cc,t = (3.7)
Il y a deux methodes numeriques pour resoudre ce probleme qui seront presentees dans les
sections 3 et 4. La premiere est la methode classique ou lequation transitoire (3.2) est resolue par
la methode des elements finis en utilisant le schema implicite de Newmark. La solution stationnaire
sera obtenue apres certains pas de temps de calcul.
La deuxieme methode est proposee pour pouvoir trouver directement la solution stationnaire
comme un probleme statique (3.6). La question est : peut on resoudre lequation stationnaire (3.6)
en utilisant la methode des elements finis ? On constate que ca marche tres bien dans le cas
subsonique. En revanche, dans le cas supersonique, cest a dire < 0, la nature de lequation (3.6)
est changee. Elle devient une equation du type hyperbolique qui nassure pas la positivite de la
matrice de rigidite et la methode des elements finis usuelle nest plus valable. Une autre equation
stationnaire qui peut etre appliquee avec la methode des elements finis sera introduite dans la
section 4. Elle est deduite a partir de la discretisation en temps de lequation transitoire par la
methode -generalisee.
On cherche la solution de lequation (3.6) sous la forme u(x) = Ae x en resolvant son equation
caracteristique : 2 v + k = 0 ce qui donne :
p p
v + v2 + 4k v v2 + 4k
1 = et 2 = (3.8)
2 2
Sur les intervalles ou u0 (x) est de signe constant, on peut representer la solution u(x) par :
Les inconnues A, B dans chaque intervalle sont determinees a partir des conditions :
(i) la continuite de u(x) et de sa derivee :
# u(x)$
# (x) =u(x)
0
x $ =0
x [, +]
x [, 0) (0, +]
(3.10)
(3.11)
62 3. Milieu a comportement uniaxial soumis a une charge mobile
A. Cas 1 : c t > 0.
Dans ce cas, cc ct > v, la vitesse de la charge est subsonique par rapport aux deux types
dondes. La solution u(x) est sous la forme de deux fonctions exponentielles qui tendent vers zero
a linfini :
f0
u(x) = exp(2c x) si x > 0
t 1t c 2c
(3.13)
f0
u(x) = exp(1t x) si x 0
t 1 c 2c
t
Demonstration.
Comme cc ct > v, v2 + 4k 0 et 1 , 2 sont reels et 1 > 0 > 2 pour tout . La condition 3.11
montre que sur un intervalle ne contenant pas 0, la derivee de u(x) ne peut changer de signe quen passant
en un point x0 ou u0 (x0 ) = 0. Si un intervalle [x0 , x1 ] est tel que u0 (x0 ) = u0 (x1 ) = 0, on a u = 0 sur cet
intervalle. Sil nexiste quun seul point x0 ]0, +[ ou x0 ] , 0[ ou u0 (x0 ) = 0, la condition dexistance
dun point x0 ou la derivee de u(x) :
sannule est AB > 0 puisque 1 2 < 0. Mais la fonction u(x) qui verifie cette condition tend vers
quand x . Donc u0 (x) ne change pas de signe sur ] , 0[ et ]0, +[ . La solution generale de u(x) se
represente par :
% &
La condition (iii) impose A = D = 0 et (i) ( u(0) =0) nous donne : B = C. La solution u(x) se simplifie
en :
% &
Il nous reste une condition (0)x u(0) = f0 pour determiner B :
On peut voir immediatement que B > 0 puisque > 0 et 1 > 0 > 2 . En consequence :
x u(0+ ) = B2 < 0 (0+ ) = c 1 (0+ ) = 1c ; 2 (0+ ) = 2c
x u(0 ) = B1 > 0 (0 ) = t 1 (0 ) = 1t ; 2 (0 ) = 2t
Le coefficient B est donc completement determine et on en deduit (3.13)
3.2 Solution analytique du probleme stationnaire 63
v 2
B. Cas 2 : 0 > c > t > 4k
Dans ce cas, v > cc ct , la vitesse de la charge est supersonique par rapport aux deux types
dondes. La solution u(x) est alors determinee comme suit :
u(x) = 0 x [0, +]
f0
u(x) = u0 (x) = p [exp(1c x) exp(2c x)] x [x0 , 0]
v2 + 4c k (3.20)
u0 h t t t t
i
u(x) = u1 (x) = t 2 exp[ 1 (x x 0 )] + 1 exp[ 2 (x x 0 )] x [, x0 ]
1 2t
ln 1c ln 2c
ou x0 = < 0 et u0 = u(x0 ).
1c 2c
Demonstration.
2
Comme 0 > c > t > 4kv , v2 + 4k 0 et 1 , 2 sont toujours reels. 2t > 1t > 0 et 2c > 1c > 0. On
utilise les memes etapes que pour le cas precedent :
Si [x0 , x1 ] (] , 0) (0, +[) tel que u0 (x0 ) = u0 (x1 ) = 0 A = B = 0 si x [x0 , x1 ] puisque
1 , 2 sont reels.
Si x0 (0, +[ tel que u0 (x0 ) = 0 A = B = 0 si x [x0 , +] car 1 , 2 > 0 et (iii) u(x) = 0
si x [x0 , +[.
la derivee u0 (x) ne change pas de signe dans (0, +[ et il y a au maximum un point x0 ] , 0) tel
que u0 (x0 ) = 0. La solution generale u(x) est donnee par :
u(x) = Ae1 x + Be2 x si x 0 (3.21)
u(x) = C 1 e 1 x + D 1 e 2 x si x [x0 , 0] (3.22)
u(x) = C 2 e 1 x + D 2 e 2 x si x [, x0 ] (3.23)
u(+) = 0 (iii) A = B = 0 car 1 , 2 > 0 u(x) = 0 si x > 0.
u(0+ ) = u(0 ) C1 = D1 .
% &
(0)x u(0) = f0 [C1 (0 )1 (0 ) C1 (0 )2 (0 )] = f0
f0
C1 = >0 puisque < 0 et 2 > 1 .
(0 )[1 (0 ) 2 (0 )]
Comme x u(0 ) = C1 (1 2 ) < 0 au voisinage de x = 0, on a (0 ) = c , 1 = 1c , 2 = 2c et on peut
en deduire :
p
v2 + 4c k f0
1c 2c = C1 = p (3.24)
c 2
v + 4c k
La fonction u(x) dans lintervalle [x0 , 0] est donc determinee par :
f0
u(x) = p [exp(1c x) exp(2c x)] (3.25)
v2 + 4c k
On rappelle que cette solution u(x) est seulement valable pour la partie au voisinage du point x = 0 et
quelle nest plus valable si sa derivee change de signe. Pour ce cas, on peut montrer facilement quil existe
ln c ln c
effectivement un point maximum qui se trouve a x0 = 1c c 2 < 0. La solution u(x) a la gauche de x0
1 2
est une fonction croissante qui sexprime par :
u(x) = C2 exp(1t x) + D2 exp(2t x) (3.26)
A partir des conditions u(x0 ) = u0 et x u(x0 ) = 0, on deduit :
u0 2t u0 1t
C2 = exp[1t x0 ] et D2 = exp[2t x0 ] (3.27)
1t 2t 1t 2t
64 3. Milieu a comportement uniaxial soumis a une charge mobile
v 2
C. Cas 3 : 4k > c > t
La vitesse de charge est supersonique par rapport aux deux types dondes. On peut expliciter
la solution generale dans ce cas comme suit :
u(x) = 0 x [0, +]
2f0
u(x) = u0 (x) = p ec x sin (c x) x [x1 , 0]
|v2 + 4c k|
(3.28)
a (xxk1 ) a
u(x) = uk (x) = uk1 e cos [a (x xk1 )] sin [a (x xk1 )]
a
x [xk , xk1 ] ; k = 2, 3, .., +
ou :
1 c 2f0
x1 = arctan ; u1 = p ec x1 sin c x1 (3.29)
c c |v2 + 4c k|
a
xk = xk1 + ; uk = uk1 e a (3.30)
a
p
v |v2 + 4a k|
a = et a =
2a 2a
avec a := t si k est pair et a := c sinon (3.31)
Demonstration.
2
Comme 4kv > c > t , v2 + 4k < 0 et 1 , 2 sont imaginaires. On va ecrire 1 et 2 sous la forme :
p
1 = + i et 2 = i avec = v /2 > 0 et = |v2 + 4k|/2.
Sur ]0, +[ , si [xi , xi+1 ] est un intervalle ou u0 (x) est de signe constant avec u0 (xi ) = u0 (xi+1 ) = 0, ces
deux points sont deux racines de u0 (x) = 0 :
On peut montrer que |u(xi+1 )| |u(xi )| comme > 0 et ui+1 , ui ne peuvent pas tendre vers zero quand
x +. Donc, A = B = 0 dans cet intervalle.
Sil y a un seul point x0 ]0, +[ ou u0 (x0 ) = 0, A = B = Cte dans [x0 , +[ u(x) quand
x + puisque > 0 u0 (x) ne change pas de signe dans ]0, +[ .
La solution generale u(x) devient :
2f0 p
u(x) = p e[(0 )x] sin [(0 )x] puisque = |v2 + 4k|/2
|v2 + 4(0 )k|
3.2 Solution analytique du probleme stationnaire 65
2f0 (0 )
x u(0 ) = < 0 comme (0 ) < 0 donc :
|v2 +4(0 )k|
2f0
u(x) = p ec x sin (c x) (3.37)
|v2 + 4c k|
On peut continuer pour trouver facilement la solution dans la zone de traction a cote en ecrivant les conditions
de continuite au point x1 (u(x
1 ) = u1 et x u(x1 ) = 0) et on obtient :
t (xx1 ) t
u(x) = u1 e cos t (x x1 ) sin t (x x1 )
t
qui va atteindre un point minimum quand : x = x2 = x1 + t . La valeur de u(x) a ce point x2 est :
u2 = u1 et (x2 x1 ) .
v 2
D. Cas 4 : 0 > c > 4k > t
Comme dans le cas precedent, la vitesse de la charge dans ce cas est supersonique par rapport
aux deux types dondes (traction et compression).
u(x) = 0 x [0, +]
f0
u(x) = u0 (x) = p [exp(1c x) exp(2c x)] x [x0 , 0]
2
v + 4c k
t (xx0 ) t (3.38)
u(x) = u1 (x) = u0 e cos [t (x x0 )] sin [t (x x0 )] x [x1 , x0 ]
t
u1
u(x) = u2 (x) = c [ c exp(1c (x x1 )) 1c exp(2c (x x1 ))] x [, x1 ]
1 2c 2
p
v |v2 + 4t k|
ou : t = ; t =
2t 2t
ln 1c ln 2c t
x0 = c c < 0 ; x 1 = x0 + ; u0 = u0 (x0 ) ; u1 = u1 (x1 )
1 2
Demonstration.
La solution de ce cas est deduite directement a partir des solutions du cas 2 et du cas 3 car les conditions
sur c et t sont identiques aux cas 2 et 3, respectivement. La solution u(x) = 0 quand x > 0 est evidente.
Si x < 0 la solution est determinee a partir du point x = 0 de droite a gauche :
x [x0 , 0], la barre est en compression, la solution u(x) = u0 (x) est la meme que celle du cas 2.
x [x1 , x0 ], on passe de la compression a la traction a partir de x = x0 , la solution u(x) = u1 (x) est
determinee sous la forme de celle du cas 3.
x [, x1 ], la barre revient en compression au point minimum x = x1 . La solution u2 (x) est
determinee de facon similaire a la solution u1 (x) sauf que londe consideree est en compression. La
derivee de cette solution ne change plus de signe jusqua donc il na plus dautres cas de solution.
66 3. Milieu a comportement uniaxial soumis a une charge mobile
v 2
E. Cas 5 : c > 0 > t > 4k
Demonstration.
La vitesse de la charge est supersonique par rapport a ct et est subsonique par rapport de cc , alors v2 +
4c,t k 0. 1 et 2 sont donc toujours reels et on peut en deduire que 2t > 1t > 0 et 1c > 0 > 2c .
c
u(+) = 0 la solution a droite doit etre en compression : u(x > 0) = Ae 2 x A > 0 comme
c
u0 (x > 0) = 2c Ae2 x < 0.
La solution a la gauche du point 0 est soit en compression soit en traction :
c
Si la solution reste en compression u(0 ) = Be1 x comme u() = 0 et comme cette solution na
c
aucun point minimum. La continuite au point x = 0 donne B = A u0 (0 ) = A1c e1 x > 0
impossible.
t t
Si la solution est en traction, elle secrit : u(x) = Ce1 x + De2 x . Les conditions de continuite nous
donnent :
A=C +D
c 2c A t (1t C + 2t D) = f0
On na que deux equations pour trois inconnues, donc le probleme a une infinite de solutions.
On peut analyser cette solution plus en detail pour montrer que le coefficient A doit verifier la condition
suivante :
f0 f0
c A (3.39)
c 2 c 2 t 2t
c
v 2
F. Cas 6 : c > 0 > 4k > t
Il ny a pas de solution
Demonstration.
c
v2 + 4c k 0 1c , 2c sont reels et 1c > 0 > 2c u(x) = Ae2 x si x > 0
v2 + 4t k 0 1t , 2t sont imaginaires.
Cette solution a toujours un point minimum x = x0 ou sa derivee change de signe (traction compression).
La solution u(x) pour x < x0 a donc forcement la forme suivante :
c c
u(x) = C2 e1 x + D2 e2 x (3.41)
Mais cette fonction ne peut pas sannuler en 2 points differents on ne peut pas avoir de solution dans ce
cas.
Les illustrations pour les solutions analytiques dans les 6 cas consideres sont presentees dans
la figure 3.3. Pour le cas 5, on trace deux solutions possibles.
3.2 Solution analytique du probleme stationnaire 67
u0
Cas 2
Cas 1
u(x)
u(x)
x0
x x
u1 u0
Cas 3 Cas 4
u3
uk
u(x)
u(x)
xk xk1 x3 x2 x1 x1 x0
uk1
u2 u1
x x
u0
Cas 6
Cas 5
A = Amin
u(x)
u(x)
x x0
Cas 5
A = Amax
?
u(x)
x u(x+
0) 6= 0 x u(x
0) =0
x
x
Remarque 3.2.1
Les cas 5 et 6 montrent que les hypotheses utilisees dans ce chapitre ne sont pas encore suffisantes pour
determiner les regimes permanents dans ces situations quand la vitesse de la charge est entre les valeurs de c c
et ct . Ces regimes permanents, sils existent, contiennent des ondes de choc qui doivent etre etudiees en tenant
compte aussi des equations en grande deformation ecrites aux surfaces de discontinuite. Pour simplifier le
probleme, on se limitera a letude des non-linearites materielles provenant des lois de comportement et on
netudiera pas les non-linearites geometriques provenant de grandes deformations.
68 3. Milieu a comportement uniaxial soumis a une charge mobile
Remarque 3.2.2
Dans le cas lineaire, Ec = Et , le probleme se simplifie puisque les cas 4, 5, 6 nexistent plus. Les solutions
u(x) dans les trois premiers cas se simplifient comme suit :
1. Si v < c
f0
u(x) = p e 1 x si x<0 (3.42)
2 + 4k
f0
u(x) = p e 2 x si x>0 (3.43)
2 + 4k
2. Si v > c et 2 + 4k 0
f0
u(x) = p (e1 x e2 x ) si x<0 (3.44)
2 + 4k
u(x) = 0 si x>0 (3.45)
3. Si v > c et 2 + 4k 0
p
2f0
2 x |v2 + 4k|
u(x) = p e sin( x) si x<0 (3.46)
|v2 + 4k| 2
u(x) = 0 si x>0 (3.47)
On va chercher la solution du probleme dynamique (3.2) par la methode des elements finis en
dynamique classique [51][115]. La solution dans le regime permanent sera obtenue en passant le
regime transitoire. On note = [L, L] le domaine de calcul en espace x et T = [0, t max ] lintervalle
de calcul en temps. Le probleme (3.2) est pose ainsi :
Chercher u(x, t) : T R tel que :
0
S u(x, t) + u(x, t) ESu0 (x, t) + ku(x, t) = f0 (x vt) (x, t) T
(3.48)
u(L, t) = u(L, t) = 0 t T
u(x, 0) = u(x, 0) = 0 x
ou on a utilise la notation u, u pour les derivees des 1 e et 2e ordre en temps et u0 pour la derivee
en espace.
Formulation variationnelle. On introduit les deux espaces des fonctions solutions et des fonc-
tions tests :
ou n designe le nombre total de noeuds, N sont les fonctions dinterpolation, u et w sont les
vecteurs des valeurs de u(x, t) et w(x) aux noeuds. On peut en deduire les approximations de la
vitesse u(x, t) et de lacceleration u(x, t) :
X
u(x, t) = NI (x)uI (t) (3.54)
I=1,n
X
u(x, t) = NI (x)uI (t) (3.55)
I=1,n
En appliquant ces approximations dans (3.51), le probleme se reduit a un systeme lineaire dequations
aux derivees partielles du 2e ordre avec les conditions initiales :
Mu + Cu + (Kb + Kr )u = F(t) (3.56)
u(0) = u(0) = 0 (3.57)
ou on appelle M la matrice de masse, C la matrice damortissement, K b , Kr les matrices de rigidite
(de la barre et des ressorts, respectivement) qui sont definies par :
Z Z
t
M= S N N dx ; C= Nt N dx (3.58)
Z Z
A. Probleme lineaire.
Remplacons (3.61)(3.62) dans (3.63), les inconnues u n+1 et un+1 seliminent et lacceleration a
linstant tn+1 peut etre exprimee en resolvant un systeme dequations lineaires :
1 n
2 t2
un+1 = M + 1 t C + (Kb + Kr ) Fn+1 + C [un + (1 1 )t un ]
2
(1 2 )t2
+(Kb + Kr ) un + t un + un (3.65)
2
70 3. Milieu a comportement uniaxial soumis a une charge mobile
Une fois que le vecteur dacceleration un+1 est determine, le deplacement et la vitesse sont calcules
directement a partir de (3.61)(3.62).
B. Probleme nonlineaire.
(k+1)
et lincrement de un+1 entre deux iterations (k, k + 1) sexprime par :
" (k)
#1
(k+1) Rn+1 (k)
un+1 = Rn+1 (3.70)
u
(k)
Il nous faut determiner la derivee de R n+1 . A partir de (3.61)(3.62), on a :
un+1 2
= (3.71)
un+1 2 t2
un+1 un+1 un+1 21
= = (3.72)
un+1 un+1 un+1 2 t
(k)
(k) Rn+1
On derive (3.68) et on obtient la matrice tangente de rigidite K t = u comme suit :
(k) 2 21 (k)
Kt = M+ C + Kbt + Kr (3.73)
2 t2 2 t
La mise a jour des solutions dans la nouvelle iteration (k + 1) est calculee par :
(k+1) (k) (k+1)
un+1 = un+1 + un+1 (3.74)
(k+1) (k) 21 (k+1)
un+1 = un+1 + un+1 (3.75)
2 t
(k+1) (k) 2 (k+1)
un+1 = un+1 + un+1 (3.76)
2 t2
3.3 Solution numerique du probleme transitoire 71
Remarque 3.3.1
1. Si on utilise des elements lineaires a deux noeuds, le vecteur dinterpolation dun element [x i , xi+1 ]
est : N =< xi+1h x , xxh
i
>. Dans le cas lineaire, les matrices elementaires Me , Ce , Keb , Ker sont
explicitees par :
Sh 2 1 h 2 1
Me = ; Ce = (3.77)
6 1 2 6 1 2
ES 1 1 kh 2 1
Keb = ; Ker = (3.78)
h 1 1 6 1 2
Dans le cas nonlineaire, Keb doit etre recalcule a chaque iteration (k) car il depend de E : E = Ec si
(k) 0 et E = Et sinon ou la deformation (k) est determinee par :
1
(k) = N(e)
,x u
e(k)
= (ui+1 ui ) (3.79)
h
2. Les parametres 1 , 2 sont choisis selon les cas de calcul. Si 1 = 2 = 0, on a un schema explicite et
si 1 = 2 = 1 le schema devient implicite total. Il est bien connu que la condition pour que le schema
de Newmark soit stable inconditionnement est 2 1 0.5. Pour des 1 , 2 quelconques, il faudra
choisir le pas de temps t suffisamment petit pour avoir la stabilite.
On va calculer quelques exemples numeriques du probleme dune barre soumise a une force
constante qui se deplace avec differentes vitesses. La barre est fixee aux deux extremites et sa
longueur varie selon les cas de calcul de telle facon quelle peut representer une barre infinie.
Les caracteristiques physiques et geometriques qui representent une structure de ballast-sol sont
donnees dans le tableau 3.1. On note que la rigidite du sol (k) est determinee par linverse du
deplacement moyen a la surface dun demi espace homogene (dont E = 100 MPa) soumis a une
force unite tangente uniforme appliquee sur une region de sa surface (dont le diametre est egal a
la largeur de la couche de ballast). La valeur de la charge utilisee dans le calcul est prise egale a
1/10 de lamplitude de la charge verticale (17 tonnes) des essieux sur la voie. On va considerer
dans un premier temps le cas lineaire et ensuite le cas nonlineaire. Dans les deux cas, la force est
deplacee a partir dune certaine position dans la barre jusqua la position ou le probleme devient
en regime permanent. Ces resultats finaux en regime permanent seront compares avec les solutions
analytiques presentees dans la section 2.
A. Probleme lineaire
Le module dYoung de la barre est constant E = 200 10 6 N.m2 . La vitesse critique corres-
pondant a ce module elastique vaut : c 333 m.s 1 . Les parametres du calcul numerique sont
donnes dans le tableau 3.2.
72 3. Milieu a comportement uniaxial soumis a une charge mobile
6
x 10
8
v = 80 m.s1 numrique
N e = 80 analytique
6
1 = 2 = 0.5
deplacement (m)
2
20 15 10 5 0 5 10 15 20
x (m)
0.5
erreur
0
0.05
erreur
0
0.5
0.05
0.1 0.15 0.2
temps (s)
1
0 0.05 0.1 0.15 0.2 0.25
temps (s)
Les resultats numeriques sont presentes dans les figures (3.4) et (3.5) en comparant avec les
courbes analytiques a la fin du calcul. Les differences entre la solution transitoire et la solution
analytique stationnaire sont aussi tracees pour voir si on est deja dans le regime stationnaire. Cette
erreur est exprimee simplement par :
P 1
N 2
|ui (t) uanalytique (xi
i=1 )|2
e= P 1 (3.80)
N
|u (x )| 2 2
i=1 analytique i
3.3 Solution numerique du probleme transitoire 73
Les pas de temps pour les calculs sont choisis tels que la charge ne peut pas se mouvoir plus loin
que la taille dun element durant un pas de temps. Pour les exemples presentes ici, on utilise un
pas de temps t qui vaut :
1 h
t = (3.81)
2 max(v, c)
5
x 10
2
numrique v = 400 ms1
1.5 analytique N e = 100
1 = 2 = 0.5
deplacement (m)
0.5
0.5
1
50 40 30 20 10 0 10 20 30 40 50
x (m)
0.4
0.2
0.2
erreur
0.02
erreur
0.4
0.01
0.6 0
0.1 0.12 0.14 0.16
0.8 temps (s)
1
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
temps (s)
B. Probleme nonlineaire
La barre se comporte avec deux modules elastiques differents : E c = 200 106 N.m2 pour
la compression et Et = 20 106 N.m2 pour la traction. Les vitesses critiques correspondant a
ces deux modules sont cc = 333.3 m.s1 et ct = 105.4 m.s1 . Les parametres de calcul utilises
sont donnes dans le tableau 3.3. Les resultats sont presentes dans les figures (3.6) et (3.7). La
difference entre les deux pentes de deplacement en compression et traction dans le cas subsonique
(3.6) montre bien leffet de la rigidite unilaterale.
Sur la figure (3.7), la courbe en pointillee represente la solution du cas lineaire pour E = E c . On
constate que les parties des solutions qui se trouvent devant la force sont identiques. La solution
qui est derriere la force dans le cas nonlineaire est plus basse frequence par rapport a celle du
cas lineaire. On constate que pour le cas supersonique, on a besoin de 280 increments de pas de
temps pour depasser le regime transitoire et dans chaque increment on doit realiser une procedure
iterative sil y a une nonlinearite (avec 2 ou 3 iterations pour notre exemple). Le cout du calcul
74 3. Milieu a comportement uniaxial soumis a une charge mobile
6
x 10
15
v = 80 ms1 numrique
N e = 160 analytique
10 1 = 2 = 1.0 lineaire
deplacement (m)
5
20 15 10 5 0 5 10 15 20
x (m)
0.5
0
erreur
0.02
erreur
0
0.5
0.02
0.1 0.15 0.2
temps (s)
1
0 0.05 0.1 0.15 0.2 0.25
temps (s)
deviendra tres cher (meme impossible !) si on doit utiliser cette procedure pour les structures en 3
dimensions.
La figure 3.8 represente le resultat dans le cas ou la vitesse de la charge a une valeur entre les
deux vitesses cc et ct . Les exemples numeriques montrent que la solution ne tend pas vers zero de
facon reguliere mais avec des chocs aux positions ou la deformation change le signe. En plus, il est
tres difficile dobtenir une stabilisation de la solution.
3.4 Solution numerique du probleme stationnaire 75
5
x 10
2
numrique v = 400 ms1
1.5 analytique N e = 200
lineaire
1 = 2 = 0.5
deplacement (m)
1
0.5
0.5
1
50 40 30 20 10 0 10 20 30 40 50
x (m)
0.2
0
0.4
erreur
erreur
0.005
0.6
0.01
0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17
0.8 temps (s)
1
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
temps (s)
v = 200ms1
N e = 200
deplacement (m)
10
50 40 30 20 10 0 10 20 30 40 50
x (m)
Considerons dabord le probleme lineaire (dont la loi de comportement secrit : = E) dune
barre de longueur 2L, lequation dynamique (3.2) secrit :
On cherche a trouver la solution stationnaire de cette equation par une methode delements finis
statique. Pour faire cela, on utilise dabord la methode directe de lintegration temporelle qui
va nous donner une equation de difference finie en temps t. La methode -generalisee est utilisee.
Ensuite, en passant au repere mobile ou les solutions stationnaires deviennent statiques, la condition
initiale est remplacee par une relation entre la solution au point considere et celle a un autre point
qui se trouve a une certaine distance. On va obtenir une equation statique qui represente lequilibre
pour chaque couple de deux points en espace x. Le developpement en detail est presente dans les
parties suivantes.
Les inconnues a linstant n + t sont determinees par des interpolations lineaires entre deux pas
de temps tn , tn+1 :
ou deux coefficients dinterpolation differents m , f sont utilises pour la force dinertie et les
autres. On utilise ensuite lapproximation de Newmark a linstant t n+1 :
t2 2 t2
un+1 = u(x, tn+1 ) un + tun + (1 2 ) un + un+1 (3.90)
2 2
un+1 = u(x, tn+1 ) un + (1 1 )tun + 1 tun+1 (3.91)
En remplacant ces discretisations (3.86)-(3.91) dans (3.85), on obtient une equation qui permet de
determiner les solutions pour chaque pas de temps t n+1 a partir de celles du pas precedent t n :
2 t2 2 t2
(1 m )S + (1 f )1 t + (1 f ) k un+1 (1 f ) ESu00n+1
2 2
2 t2
(1 m )S + (1 f )1 t f k un
2
t (3.92)
(1 m )S [2 2(1 f )1 ] t un
2
t2 2 t2
[(1 2 m )S (1 f )(2 1 )t ] un f ESu00n
2 2
2 t2
= f0 [(1 f )(x vtn+1 ) + f (x vtn )]
2
Demonstration.
On va eliminer les variables en tn+t dans (3.85) en utilisant (3.86)-(3.89) et puis un+1 et un+1 pour
1
la methode -generalisee (generalized-) est proposee par Chung et Hulbert [22] a partir de la methode de
Hilber-Hughes-Taylor (HHT-) [50]
3.4 Solution numerique du probleme stationnaire 77
obtenir une equation dont les variables ne sont que un+1 et un , un , un . A partir de (3.90)(3.91), on peut
exprimer un+1 , un+1 en fonction de un+1 et des variables a linstant tn :
2 2 2 1 2
un+1 = un+1 un un un (3.93)
2 t2 2 t2 2 t 2
21 21 2 21 (2 1 )t
un+1 = un+1 un + un + un (3.94)
2 t 2 t 2 2
On remplace un+1 , un+1 dans lequation dont les variables en tn+t sont deja eliminees :
S[(1 m )un+1 + m un ] + [(1 f )un+1 + f un ] ES[(1 f )u00n+1 + f u00n ]
+k[(1 f )un+1 + f un ] = f0 [(1 f )(x vtn+1 ) + f (x vtn )] (3.95)
et on obtient lequation suivante apres quelques simplifications :
2 21
(1 m ) S + (1 f ) + (1 f )k un+1 (1 f )ESu00n+1
2 t2 2 t
2 21
(1 m ) S + (1 f ) f k un
2 t2 2 t
2 21
(1 m ) S 1 (1 f ) un
2 t 2
1 m 2 1
1 S (1 f ) t un ESf u00n
2 2
= f0 [(1 f )(x vtn+1 ) + f (x vtn )] (3.96)
2 t2
En multipliant cette equation par 2 , on obtient (3.92).
Dans le regime permanent, on a toujours une relation supplementaire entre les solutions aux
deux pas de temps (tn ,tn+1 ) :
Passons au repere mobile defini par x = xvtn+1 , les vitesses et accelerations dans (3.92) peuvent
sexprimer par :
un un+1 (x + vt)
un = v = v (3.98)
x x
2
un 2
2 un+1 (x + vt)
un = v2 = v (3.99)
x2 x2
Lequation (3.92) devient stationnaire :
2 t2 2 t2
(1 m )S + (1 f )1 t + (1 f ) k u(x) (1 f ) ESu00 (x)
2 2
2 t2
(1 m )S + (1 f )1 t f k u(x + vt)
2
t
+ (1 m )S [2 2(1 f )1 ] vt u0 (x + vt)
2
v 2 t2 00
[(1 2 m )S (1 f )(2 1 )t ] u (x + vt)
2
2 t2
f ESu00 (x + vt)
2
2 t2
= f0 [(1 f )(x) + f (x + vt)] (3.100)
2
78 3. Milieu a comportement uniaxial soumis a une charge mobile
On voit que si v = 0, cette equation se reduit a lequation statique : ku(x) ESu 00 (x) = f0 (x).
Afin de reduire lecriture dans le calcul, on reecrit (3.100) sous la forme plus simple suivante :
ou :
2 t2
A1 = (1 m )S + (1 f )1 t + (1 f ) k (3.102)
2
2 t2
A2 = (1 f ) ES (3.103)
2
2 t2
B1 = (1 m )S + (1 f )1 t f k (3.104)
2
t
B2 = (1 m )S [2 2(1 f )1 ] vt (3.105)
2
v 2 t2 2 t2
B3 = [(1 2 m )S (1 f )(2 1 )t ] + f ES (3.106)
2 2
2 t2
F1 = (1 f ) f0 (3.107)
2
2 t2
F2 = f f0 (3.108)
2
A. Formulation variationnelle
Comme il y a une quantite vt qui intervient dans cette equation, on va definir un autre espace
qui est defini par : = [L, L + vt]. On introduit deux espaces : U lespace des fonctions
solutions et W lespace des fonctions tests qui sont definis par :
U = {u | u H 1 (), u(L) = u(L + vt) = 0} (3.110)
1
W = {w | w H (), w(L) = w(L) = 0} (3.111)
La formulation variationnelle de (3.109) secrit comme suit :
Chercher u(x) U : R tel que w W :
(A u(x), w(x)) + (A u0 (x), w0 (x)) (B u(x + vt), w(x)) + (B u0 (x + vt), w(x))
1 2 1 2
0 0
+ (B3 u (x + vt), w (x)) = (F1 (x) + F2 (x + vt), w(x))
(3.112)
R
ou la notation (, ) designe : (, ) = (, ) dx. Les conditions aux limites sont incluses dans cette
formulation en utilisant le fait que :
u00 (x), w(x) = u0 (x), w(x) |L 0 0
L u (x), w (x)
= u0 (x), w0 (x) (3.113)
3.4 Solution numerique du probleme stationnaire 79
S
Le domaine peut ensuite etre discretise par un ensemble de sous domaines : = e . Dans
e
chaque sous domaine e , u et w sont determines en fonction de leurs valeurs aux noeuds grace aux
fonctions dinterpolation :
ou ne est le nombre de noeuds de lelement e . Ici, on utilise la methode de Galerkin standard qui
suppose que N N.
On discretise le domaine [L, L] avec N e elements. Pour chaque element e defini dans [x i , xi+1 ],
on a :
Z xi+1 Z xi+1
h h e
a(u , w ) = wJ [ A1 NJ (x)NI (x)dx + A2 N0J (x)N0I (x)dx
xi xi
| {z } | {z }
Mea Kea
Z xi+1 Z xi+1
B1 NJ (x)NI (x + vt)dx + B2 NJ (x)N0I (x + vt)dx
xi x
| {z } | i {z }
Meb Ceb
Z xi+1
+ B3 N0J (x)N0I (x + vt)dx]uI
xi
| {z }
Keb
= wJ KeJI uI (3.117)
Z xi+1
f (wh )e = wJ [ NJ (x)(F1 (x) + F2 (x + vt))dx]
xi
| {z }
Fe
= wJ FeJ (3.118)
Ku = F (3.120)
La procedure delements finis peut etre realisee normalement avec les matrices elementaires
definies dans (3.117). Les matrices M ea , Kea sont calculees sans difficulte comme dans le cas classique.
En revanche, les integrations de M eb , Ceb , Keb sont anormales comme la fonction u(x + vt) a un
decalage en variable par rapport a w(x). Ces integrations ont les memes formes et peuvent etre
calculees de la meme facon. On evalue par exemple une integration sur un element e qui se trouve
dans [xi , xi+1 ] comme suit :
Z xi+1
I= w(x)u(x + vt)dx (3.121)
xi
On rappelle lhypothese qui impose que vt soit toujours inferieure a la taille dun element (vt <
xi+2 xi+1 ) et on peut en deduire :
(3.124)
u(x + vt)
1
w(x) 2
1 2
Si on discretise la barre par des elements dont les tailles sont regulieres h et si dans chaque
element ([xi , xi+1 ]), on utilise une fonction dinterpolation lineaire :
xi+1 x x xi
N(x) = < , > (3.125)
h h
on peut determiner explicitement les matrices elementaires M ea , Kea , Meb , Ceb , Keb , Fe dun element
e comme suit :
e A1 h 2 1 0 e A2 1 1 0
Ma = ; Ka = (3.126)
6 1 2 0 h 1 1 0
3.4 Solution numerique du probleme stationnaire 81
B1 h (1 )2 (2 + ) (1 + 3 2 3 ) 3
Meb = 3 2 2 (3.127)
6 (1 ) (2 )(1 + 2 2 ) (3 )
B2 (1 2 ) (1 2 2 ) 2
Ceb = 2 (3.128)
2 (1 ) 1 + 2( 2) ( 2)
B3 1 2 1
Keb = (3.129)
h 1 1 2
e=ef F2 e6=ef 0
F = ; F = (3.130)
F1 + (1 ) F2 0
vt
ou = h 1 et ef est lindice de lelement qui se trouve juste devant la position de la force.
Une equation par elements finis est toujours equivalente avec un certain schema de difference
finie. Dans ce paragraphe, on va essayer decrire lequation en difference finie equivalente a (3.120)
pour minimiser les erreurs numeriques causees par les approximations en reglant les parametres
m , f , 1 , 2 . Comme les matrices elementaires secrivent sur 3 noeuds, le schema en difference
finie equivalent peut etre ecrit sous la forme :
aui1 + bui + cui+1 + dui+2 = 0 (3.131)
On dit que ce schema est convergent si le developpement en serie de la fonction :
f (h) = a(h)u(x h) + b(h)u(x) + c(h)u(x + h) + d(h)u(x + 2h) (3.132)
tend vers zero en utilisant lequation stationnaire quand h 0. On va simplifier les expressions de
a, b, c, d en prenant un pas de temps t qui vaut : t = hv (en consequence = 1) et qui donne
la matrice de rigidite elementaire suivante :
A h A A1 h A2 B1 h B2 B3 B1 h B2 B3
1 2
+ + +
3 h 6 h 3 2 h 6 2 h
Ke = (3.133)
A1 h A2 A1 h A2 B1 h B2 B3 B1 h B2 B3
+ + +
6 h 3 h 6 2 h 3 2 h
On peut en deduire les coefficients a, b, c, d :
A1 h A2
a = Ke21 = (3.134)
6 h
2A 1h 2A2 B1 h B2 B3
b = Ke11 + Ke22 = + (3.135)
3 h 6 2 h
e e A1 h A2 2B1 h 2B3
c = K12 + K23 = + (3.136)
6 h 3 h
B 1 h B 2 B 3
d = Ke23 = + (3.137)
6 2 h
Le polynome (3.132) est developpe jusquau 6 e ordre en serie de Taylor avec differents parametres
m , f , 1 , 2 . Par exemple, avec les parametres 1 = 0.5, 2 = 0.5, m = 0.25, f = 0.5, le
developpement de f (h) nous donne :
h3 h4 (1) h5 (2) h5
f (h) g(x) + D g(x) + D g(x) + ku00 (x) + O(h6 ) (3.138)
4v 2 8v 2 12v 2 48v 2
h3 h h2 (2) h2 00
= 2
g(x) + D (1) g(x) + D g(x) + ku (x) + O(h ) 3
(3.139)
4v 2 3 12
82 3. Milieu a comportement uniaxial soumis a une charge mobile
ou g(x) = (E v 2 )Su00 (x) + vu0 (x) + ku(x) est la fonction de lequation stationnaire. Comme
lequation stationnaire est definie par g(x) = 0, lerreur de cette approximation (f (h) = 0) sex-
prime par :
kh2 00
e(h) u (x) + O(h3 ) (3.140)
12
Cest la meilleure approximation quon a trouvee dont lerreur est du 2 e ordre par rapport a la
taille des elements h 2 .
Si on ajoute dans f (h) une correction numerique qui est determinee par :
kh3
(ui1 2ui + ui+1 ) (3.141)
48v 2
kh5
dont le developement en fonction de h vaut , on peut supprimer lerreur dans (3.138).
48v 2
Revenons a la methode des elements finis, on peut ajouter dans la matrice de rigidite elementaire
une matrice de correction suivante pour ameliorer la precision :
kh3 1 1 0
Kcorr = (3.142)
48v 2 1 1 0
La procedure presentee ci-dessus peut etre etendue facilement aux cas nonlineaires, surtout
pour le comportement du type unilateral quon considere ici (3.3)(3.4). Dans ce cas, les matrices
tangentes elementaires peuvent etre determinees explicitement par les memes formules que (3.117).
Par exemple pour un element e, on a :
ou les trois matrices Mea , Meb , Ceb sont toujours constantes. Deux matrices K ea , Keb sont calculees
en fonction de letat de la deformation dans lelement considere.
Si on utilise des elements lineaires a deux noeuds, la deformation dans un element (e) est
constante : e = h1 (ui+1 ui ) et en consequence les coefficients A 2 , B3 sont aussi constants dans
chaque element. Les matrices Kea et Keb sont determinees par :
(e)
A 1 1 0
Kea = 2 (3.144)
h 1 1 0
(e) (e+1)
B 1 1 0 B 0
Keb = 3 + 3 (3.145)
h 1 1 0 h 0
(e) (e)
ou A2 , B3 sont determines comme dans (3.103) et (3.106) :
(e) 2 t2
A2 = (1 f ) ES (3.146)
2
(e) v 2 t2 2 t2
B3 = [(1 2 m )S (1 f )(2 1 )t ] + f ES (3.147)
2 2
2
Un exemple pour lautre choix de parametres , si 1 = 0.5, 2 = 0.5, m = 0., f = 0. (qui correspond a une
discretisation en temps par le schema Newmark), lerreur e(h) vaut :
h h2 h2 h2
e(h) S u000 (x) ES u00 (x) + v u000 (x) + ku00 (x) + O(h3 )
3 6 12 4
3.5 Validation numerique 83
On trate le meme probleme que celui du calcul transitoire qui est pose dans la section 1.2 pour
valider le schema numerique propose. Les parametres physiques et geometriques de la structure sont
donnes dans le tableau 3.1. Comme avec la validation de la methode du calcul transitoire, les deux
regimes subsoniques et supersoniques sont consideres pour les deux cas : lineaire et nonlineaire.
Dans tous les exemples suivants, on utilise les memes parametres numeriques : m = 0.25,
f = 0.5, 1 = 0.5, 2 = 0.5, = 1 (t = hv ).
On calcule dabord le cas lineaire. Les parametres sont donnes dans le tableau 3.4. Comme la
force est fixe dans ces calculs, on peut donc utiliser des longueurs de barre qui sont respectivement
les moities de celles des calculs transitoires en gardant les memes tailles des elements.
6
x 10
8
v = 80 ms1 numerique
N e = 40 analytique
6
deplacement (m)
0
10 8 6 4 2 0 2 4 6 8 10
x (m)
5
x 10
2
v = 400 ms1 numerique
e analytique
1.5 N = 100
deplacement (m)
0.5
0.5
1
25 20 15 10 5 0 5 10 15 20 25
x (m)
10
9 8 7 6 5 4 3 2 1 0 1
x (m)
Les figures (3.10) et (3.11) representent les comparaisons entre les resultats numeriques et
les solutions analytiques dans deux cas de vitesse : le cas subsonique (v = 80 m.s 1 ) et le cas
supersonique (v = 400 m.s1 ).
On constate que les resultats obtenus sont tres proches des solutions analytiques. Les erreurs
dans les cas subsoniques et supersoniques sont respectivement e = 1.2% et e = 1.8%.
On profite aussi de cet exemple pour montrer leffet de la correction numerique. Dans la figure
3.12, on trace la solution analytique en comparant le resultat numerique avec ou sans la correction.
On voit bien que la correction peut donner une amelioration importante a la precision du calcul.
Lerreur dans le calcul sans correction est e = 6.7%.
Les parametres de calcul sont donnes dans le tableau 3.5 pour les deux cas subsonique et
supersonique. Les courbes qui sont tracees dans les figures 3.13 et 3.14 montrent des resultats
assez precis par rapport aux solutions analytiques. La convergence dans le cas subsonique est
obtenue apres 4 iterations. Pour le cas supersonique, le nombre diterations necessaires est 12. On
voit bien que cette methode permet une diminution importante du volume de calcul par rapport a
la methode transitoire.
Remarque 3.5.1
Dans le cas transonique ct < v < cc , on ne peut pas trouver la solution car la procedure iterative ne converge
3.6 Conclusion 85
6
x 10
15
v = 80 ms1 numerique
N e = 100 analytique
10
deplacement (m)
5
10 8 6 4 2 0 2 4 6 8 10
x (m)
0.5
0.5
1
50 40 30 20 10 0 10 20 30 40 50
x (m)
pas puisque il ny a pas de solution en petite deformation comme on a montre dans le paragraphe precedent.
3.6 Conclusion
On a etudie dans cette partie les phenomenes physiques de la propagation dondes dans une
barre non-lineaire de type unilateral soumise a une charge mobile constante. Les solutions analy-
tiques dans les differents cas de vitesse de charge montrent les differents mecanismes de propagation
des ondes dans la barre. Elles montrent aussi que, dans le cadre des hypotheses quon a utilisees,
on ne peut atteindre que le regime permanent si la vitesse de la charge est soit inferieure a la plus
petite vitesse donde dans la barre (absolument subsonique), soit superieure a la plus grande vitesse
donde dans la barre (absolument supersonique). Autrement, letat des propagations dondes dans
la barre devient beaucoup plus complexe avec des ondes de choc qui apparaissent.
86 3. Milieu a comportement uniaxial soumis a une charge mobile
Dans ce chapitre, on a presente aussi des methodes numeriques pour resoudre le probleme
non-lineaire dune barre soumise a une charge mobile avec une vitesse constante. Les solutions
analytiques de ce probleme sont utilisees pour evaluer la precision des schemas numeriques. Dans un
premier temps, on a utilise un schema dintegration numerique classique (Newmark) pour calculer le
probleme dans le regime transitoire. Neanmoins, cette methode demande un gros maillage et aussi
un grand nombre de pas de temps de calcul si on veut acceder au regime stationnaire. En raison
de cela, une autre methode a ete proposee dans un deuxieme temps afin devaluer directement
la solution stationnaire dans un nouveau repere mobile qui sattache a la position actuelle de la
force. Au lieu de resoudre lequation dans le repere mobile qui devient une equation hyperbolique
dans le cas supersonique, on a resolu une equation elliptique equivalente a celle obtenue a partir
de lequation transitoire discretisee en temps. Les comparaisons montrent que les resultats obtenus
par ce schema saccordent bien avec les resultats analytiques dans les deux cas subsonique et
supersonique.
En conclusion :
4.1 Introduction
Dans la modelisation des voies ferrees, le rail est considere comme une poutre en flexion couplee
avec la structure dassise en dessous. Pour simplifier la complexite des structures dassise, on
suppose que chaque composante (les blochets, le ballast, les differentes couches du sol ...) est
uniforme et continue dans la direction horizontale. Cela conduit a un probleme dune poutre posee
sur un milieu multi-couche. Dans un modele plus simple, ce milieu multicouche peut etre modelise
par un systeme discret dans la direction verticale ou on definit a chaque interface un degre de
liberte (figure 4.1). Par exemple, chaque couche peut etre modelisee par un systeme simple de
masse-ressort-amortisseur. Donc, sil existe n couches differentes, on a un systeme discret de n
degres de liberte. La structure globale devient un systeme discret de n degres de liberte (poutre +
systeme discret).
z f0 (t)(x vt)
x +
wr
w1 structure
(1) dassise
w2
(2)
: systeme
: wn discret de
(n) n DDL
On cherche les reponses dynamiques de ce systeme dues a une charge ponctuelle (dont lampli-
tude f0 (t) est soit constante soit variable en fonction du temps) qui se deplace avec une vitesse v
constante. Les equations dynamiques de ce systeme peuvent etre etablies en ecrivant les equilibres
entre le rail et la structure dassise. Appelons f r la force dinteraction entre ces deux parties,
lequation dynamique du rail qui est considere comme une poutre dEuler-Bernouilli secrit :
ou r Sr ,Er Ir designent la masse et la rigidite en flexion sur une unite de longueur de la poutre, w r
est le deplacement vertical. Lacceleration de gravite g est ajoutee pour tenir compte du poids de
la poutre et des differentes parties de la structure dassise.
Le systeme discret de n degres de liberte peut toujours etre decrit par 3 matrices de dimensions
n n qui representent la masse ma , lamortissement a et la rigidite ka . En notant w le vecteur
de deplacement de ces n degres de liberte, lequation dynamique sous la sollicitation de la force f r
doit verifier le systeme dequations qui secrit sous la forme suivante :
(x, t) R R+ :
mw(x, t) + w(x, t) + kw(x, t) + jw (4) (x, t) = f 0 (t)(x vt) + g (4.4)
w(x, 0) = w 0 ; w(x, 0) = v 0
Alors, on est conduit a resoudre le systeme non-lineaire (4.4). La methode dintegration directe
ou la position de la charge change apres chaque pas de temps est toujours utilisable mais avec
un grand volume de calcul. On va presenter dans la partie suivante la methode numerique pour
resoudre le probleme dans le repere mobile. Ensuite, la solution analytique pour le cas plus simple
(4.5) sera presentee. A la fin, les resultats analytiques et numeriques sont compares pour valider le
schema numerique propose.
4.2.1 Methodologie
i.e. la fonction h(x, t) peut etre decomposee comme un produit de deux fonctions h v et ht qui ne
dependent que de la position de la charge et du temps, respectivement.
En utilisant cette decomposition, lequation stationnaire discretisee en temps sera ecrite pour
le premier terme hv (x vt) et on va obtenir un systeme dequations dont le vecteur force est fixe
en espace dans le repere mobile et qui, cette fois, est dynamique.
Comme la fonction w(x, t) dans (4.4) est un vecteur de dimension n, lapplication de (4.6) pour
w(x, t) secrit : w i (x, t) = wvi (x, t)w ti (x, t), i = 1 n. Afin de simplifier les ecritures des equations
dans le calcul qui suit (et qui doit normalement secrire avec les indices), on propose dutiliser un
operateur simple qui est defini comme suit :
c = a b ci = ai bi
(pas de somme en i) (4.7)
B = A a B ji = Aji ai
ab=ba (4.8)
A(a b) = (A a)b = (A b)a (4.9)
w(x, t) = wv wt + wv wt (4.11)
w(x, t) = wv wt + 2wv wt + wv wt (4.12)
En remplacant les 3 equations ci-dessus dans (4.4), on obtient un systeme dequations differentielles
en temps pour la variable w v :
mwv + wv + kw v + jw(4)
v = f 0 (t)(x vt) + g (4.13)
m = m wt (4.14)
= 2 m wt + wt (4.15)
k = m wt + wt + k wt (4.16)
j = j wt (4.17)
Demonstration.
On applique (4.10) dans (4.4) :
Equation discretisee en temps pour w v (x vt). De meme facon que dans le cas uniaxial,
on discretise le vecteur w v (x vt) dans lequation (4.13) par la methode -generalisee [22] et on
obtient un syteme dequations qui relie la solution a linstant t n+1 en fonction de celle a linstant
tn :
2 t2 2 t2
(1 m )m + (1 f )1 t + (1 f ) k wv(n+1) + (1 f ) jw(4)
v(n+1)
2 2
2 t2
(1 m )m + (1 f )1 t f k wv(n)
2
t
(1 m )m [2 2(1 f )1 ] t wv(n)
2
t2 2 t2
[(1 2 m )m (1 f )(2 1 )t ] wv(n) + f jw(4)
v(n)
2 2
2 t2 2 t2 2 t2
= (1 f )f (t)(x vtn+1 ) + f f(t)(x vtn ) + g
2 2 2
(4.18)
ou t = tn+1 tn . Pour simplifier les ecritures, on reecrit cette equation sous la forme :
A1 = A1 wt + A3 wt + A4 wt (4.34)
A2 = A2 wt (4.35)
B 1 = B 1 wt + B 5 wt B 8 wt (4.36)
B 2 = B 2 wt B 6 wt (4.37)
B 3 = B 3 wt B 7 wt (4.38)
B 4 = B 4 wt (4.39)
En utilisant (A b)a = A(a b), on arrive a la derniere etape ou on trouve la formule ecrite dans
le repere mobile qui est equivalente a (4.4) :
(x, t) R R+ :
A1 w(x, t) + A2 w(4) (x, t)
B 1 w(x+ , t) + v B 2 w0 (x+ , t) v 2 B 3 w00 (x+ , t) + B 4 w(4) (x+ , t)
(4.47)
+ A3 w(x, t) B 5 w(x+ , t) v B 6 w0 (x+ , t) + v 2 B 7 w00 (x+ , t)
+ [A4 w(x, t) + B 8 w(x+ , t)]
= F 1 (t) (x) + F 2 (t) (x+ ) + G
92 4. Milieu unidimensionnel en flexion soumis a une charge mobile
A ce stade, on peut ecrire la formulation variationnelle pour (4.47) et puis etablir le systeme
dequations en elements finis. Supposons que le domaine de calcul est = [L, L] avec les condi-
tions aux limites : w(x < L) = w(x > L) = w, x (x < L) = w,x (x > L) = 0, la formulation
variationnelle du probleme (4.47) secrit :
Remarque 4.2.1
1. On utilise ici une matrice N qui est construite a partir de la fonction dinterpolation de Hermite pour
la poutre (wr ) et des fonctions dinterpolation lineaires pour wj et wj+1 . Au point x [xi , xi+1 ], cette
fonction N permet devaluer le vecteur w(x, t) a partir des valeurs aux 2 noeuds (i, i + 1) comme suit :
wi (t)
w(x, t) = N(x) (4.54)
wi+1 (t)
4.2 Systeme dequations dynamiques dans le repere mobile 93
0
ou wi = {wr , r , w2 , .., wn }ti (ri = wri est la rotation de la poutre au noeud i), N(x) secrit donc sous la
forme de la matrice n 2(n + 1) suivante :
N1 N2 0 . 0 N 3 N4 . 0 0
0 0 N1j . 0 0 0 N2j . 0
N(x) = .
(4.55)
. . . . . . . . .
n n
0 0 0 . N1 0 0 0 . N2 n2(n+1)
ou :
(x xi+1 )2 [h + 2(xi x)] (x xi )(x xi+1 )2
N1 = ; N2 = (4.56)
h3 h2
(x xi )2 [h + 2(xi+1 x)] (x xi )2 (x xi+1 )
N3 = ; N4 = (4.57)
h3 h2
x i+1 x x x i
N1j = ; N2j = (4.58)
h h
2. Dans les formulations des matrices elementaires, il y a des integrations avec un decalage qui portent
sur N+ . De meme facon que pour la methode presentee dans le chapitre precedent, ces integrations
sont calculees sur deux elements et les matrices elementaires quon doit utiliser ont la dimension de
2(n + 1) 2(n + 1). Dans le cas non-lineaire, elles sont calculees par la technique dintegration des points
de Gauss.
Remarque 4.2.2
1. Le systeme dequations dynamiques dans le repere mobile peut etre resolu en utilisant un certain schema
dintegration directe (par exemple Newmark) et sa solution stationnaire est obtenue apres le regime
transitoire. On remarque que cette solution transitoire nest pas la vraie solution transitoire du probleme
(4.4). En fait, la solution statique ici est celle du probleme ou la charge est constante. La solution
transitoire equivalente a celle du probleme ou on a la structure en regime permanent due a une charge
constante et puis la charge commence a varier en fonction du temps.
2. Dans le cas ou on a une charge constante qui se deplace avec une vitesse v, lequation dans le repere mobile
est statique, la solution stationnaire peut etre trouvee directement en resolvant le systeme dequations
statique :
Kw = F (4.59)
On va etudier la precision des equations discretisees en elements finis donnees par les equations
proposees. Pour simplifier le probleme, on va considerer un probleme ou on a seulement la poutre
couplee avec les ressorts de Winkler. Lequation dynamique devient scalaire et est definie par :
Le systeme dequations des elements finis dont les matrices elementaires sont presentees dans
lannexe C nous donne a chaque noeud une equation de type difference finie comme suit :
X
am m c c k k
` wi+` + b` i+` + a` wi+` + b` i+` + a` wi+` + b` i+` = 0 (4.61)
`=1,2
ou am ,bm ,ac ,bc ,ak ,bk sont des coefficients constants qui sont evalues en fonction de ,,k,E,S et h
(h = xi+1 xi ). On dit que ce schema de difference finie est convergent si :
X h
f (h) = am m c c
` (h)w(x + `h) + b` (h)(x + `h) + a` (h)w(x + `h) + b` (h)(x + `h)
`=1,2
i
+ak` (h)w(x + `h) + bk` (h)(x + `h) 0 quand h 0 (4.62)
94 4. Milieu unidimensionnel en flexion soumis a une charge mobile
On developpe ce polynome en serie de Taylor pour mesurer lerreur pour differentes valeurs des
parametres m , f , 1 , 2 et t. Par exemple, on peut expliciter le developpement jusquau 6 e
ordre dans le cas ou m = 0.25, f = 0.5, 1 = 0.5, 2 = 0.5, t = h/v (les memes que ceux
utilises pour le cas uniaxial) avec laide de Mathematica :
On a teste avec plusieurs valeurs de m , f , 1 , 2 et on a verifie que lerreur est toujours plus
petite avec ces valeurs des parametres : m = 0.25, f = 0.5, 1 = 0.5, 2 = 0.5, t = h/v. On va
les utiliser dans tous les calculs suivants.
4.3 Validation
La validation est realisee dans un cas simple ou on considere un probleme lineaire dune poutre
posee sur la fondation de Winkler. La charge est soit constante soit harmonique. Les resultats
numeriques obtenus par la methode proposee sont compares avec ceux obtenus par une methode
semi-analytique.
On presente ici les phenomenes vibratoires dans un cas plus simple dune poutre posee sur la
fondation de Winkler soumise a une charge harmonique. La solution est obtenue par un calcul
semi-analytique. Le regime stationnaire est cherche. En notant S, EI, , k respectivement la
masse, la rigidite en flexion de la poutre, lamortissement et la rigidite du systeme des ressorts,
lequation dequilibre dynamique du probleme se represente par :
Par un changement de variable x := x vt w(x, t) := w(x vt, t), lequation dans le repere
mobile devient :
Comme la force est harmonique de frequence , la solution stationnaire de cette equation peut etre
ecrite sous la forme :
Sv 2 2 (2S + i)v k i S 2
4 + =0 (4.69)
EI EI EI
z
T (x + dx)
x
M (x) M (x + dx)
T (x) dx
Remarque 4.3.1
Afin danalyser le comportement des j de facon plus generale, on peut convertir (4.69) en une forme
r
4 4k
non-dimensionnelle en posant = avec = :
EI
1 2 2i
4 2 2 (i + ) + =0 (4.76)
4
ou on a introduit trois parametres non-dimensionnels : la frequence , la vitesse et lamortissement :
v
= , = , = (4.77)
0 v0 0
s s
k 4EIk p
0 = , v0 = 4 , 0 = 2 Sk (4.78)
S (S)2
Les signes des parties reelles et imaginaires de quatre racines ne dependent que des valeurs de , ,
. Dans le cas ou = 0 (charge constante), on peut montrer que lon peut toujours distinguer 2 ondes a
gauche et 2 autres a droite de la charge ([39], chap. 13). Si 6= 0, le probleme devient plus complique. Une
discussion detaillee peut etre trouvee dans [7].
Afin de reduire la taille du domaine de calcul mais sans perturber le resultat a cause des ondes
reflechies par les frontieres, il est necessaire dintroduire des conditions aux limites speciales aux
frontieres.
On propose ici dajouter une couche delements dont lamortissement permet dabsorber les
ondes incidentes. Les ondes qui entrent dans cette couche seront attenuees beaucoup plus rapide-
ment mais de telle facon quil ny a pas de reflexion. Notons la taille de la couche absorbante par
, on definit lamortissement par une fonction de x :
si |x| < L
(x) = (4.79)
|x| L 4
+
si |x| > L
xL 4
= +
x
L
couche absorbante
En resolvant lequation (4.66) avec une telle fonction damortissement (x), on peut verifier
que la solution na pas de discontinuite de w et de ses derivees au point x = L.
4.3 Validation 97
La figure 4.4 presente un exemple pour montrer leffet de la couche absorbante pour un probleme
ou la charge est fixe en espace. On considere une poutre (EI = 10 8 N m2 ) qui est posee sur un
systeme de ressorts uniformes (k = 10 3 N m2 ) sans amortissement. La force harmonique est ap-
pliquee sur le point x = 0 : f (t) = f0 eit , = 5.8 rad/s, f0 = 1. Un domaine qui se trouve
dans (250m, 250m) est maille par 100 elements identiques. La taille de la couche absorbante est
= 150m de chaque cote. On peut constater que les couches absorbantes attenuent les ondes qui
se propagent loin de la force et donc la solution dans la partie interessante est bien en accord avec
celle analytique. Comme lamortissement dans ce cas est nul, le resultat numerique obtenu sans
couche absorbante devient tres mauvais (la courbe pointillee dans la figure 4.4).
5
x 10
2
analytique
FEM
1.5 solution sans
condition absorbante
1
deplacement vertical (m)
0.5
0.5
1.5
couche absorbante couche absorbante
2
250 200 150 100 50 0 50 100 150 200 250
X (m)
Figure 4.4. Effet de la couche absorbante : cas dune charge fixe harmonique
Pour les exemples numeriques, on considere un rail de type UIC60 pose sur un sol solide. Les
caracteristiques du rail et du sol sont donnees dans le tableau 4.1. On determine aussi les valeurs des
parametres adimentionaux 0 , v0 , 0 (4.78) qui correspondent aux caracteristiques de ce systeme.
On presente dans la suite les calculs de validation qui sont effectues avec deux types differents
de charge : soit avec une force constante soit avec une force harmonique.
Force constante
Le regime permanent avec une charge mobile dont lamplitude est constante peut etre obtenu
par un calcul statique. La figure (4.5) presente les resultats de deplacement vertical due a une
charge de 10 kN dans deux cas de vitesse : le cas subsonique (v = 0.5 v 0 ) et le cas supersonique
(v = 1.5 v0 ).
Dans le cas subsonique, la solution est attenuee rapidement, on na besoin que de 50 elements
sur une longueur de calcul de 20 m. La couche absorbante nest pas necessaire dans ce cas. Quand
la vitesse v depasse v0 , la solution devient de type sinusodale, en plus, elle est vraiment en haute
98 4. Milieu unidimensionnel en flexion soumis a une charge mobile
frequence dans la zone devant la force. On utilise alors un maillage de 200 elements sur une longueur
de 50 m ou il y a de chaque cote une couche absorbante de longueur = 25 m. La courbe pointillee
represente le resultat sans utiliser la couche absorbante. Pour obtenir un resultat de meme qualite
avec une condition aux limites normale, on a besoin dun maillage de 500 elements sur un longueur
de 100 m.
4 v = 0.5 v 4 v = 1.5 v
cr cr
x 10 x 10
1 3
numerique numerique
analytique analytique
2
0
1
deplacement vertical (m)
1
0
1
2
3
3
4 4
10 5 0 5 10 30 20 10 0 10 20 30
X (m) X (m)
Figure 4.5. Validation du modele poutre - fondation de Winkler : cas dune force constante
Force harmonique
La validation sur le probleme dynamique est realisee avec une force harmonique dont la
frequence = 0.50 . Deux cas de vitesse (subsonique et supersonique) sont toujours etudies.
La figure 4.6 presente le resultat stationnaire pour v = 0.5v 0 et le deplacement au point x = 0
en fonction du temps. La longueur de la poutre calculee est 30 m, le nombre delements est 90. Le
pas de temps dt est pris egal a 0.05 2/. On peut avoir une solution numerique correcte sans
la couche absorbante.
Pour le cas supersonique, la solution stationnaire est obtenue plus difficilement. On a besoin
dun maillage de 500 elements dans un domaine de 100 m pour obtenir le resultat presente dans la
4.4 Conclusion 99
figure 4.7. Le pas de temps utilise est 0.02 2/ et il faut calculer 400 pas pour arriver au regime
stationnaire. On voit que les ondes se propagent par groupe et on observe leffet Doppler a cause
de la vitesse de la charge. Les epaisseurs des couches absorbantes prises pour ce calcul sont 25 m.
4 v = 0.5 v 3
cr
x 10 x 10
2 1
numerique numerique
analytique analytique
1
1
0
2
5 1
15 10 5 0 5 10 15 0 0.05 0.1 0.15
X (m) temps (s)
4 v = 1.5 vcr 4
x 10 x 10
4 1
numerique numerique
analytique analytique
3
deplacement au point x = 0 (m)
2
deplacement vertical (m)
0 0
4 1
50 0 50 0 0.1 0.2 0.3 0.4
X (m) temps (s)
4.4 Conclusion
5.1 Introduction
Le schema de ce modele est presente dans la figure 5.1. La structure est caracterisee par la masse
de section S et la rigidite en flexion EI de la poutre et aussi par la rigidite k et lamortissement
de la fondation.
z
f0 (t)(x vt)
x +
(x, t) R R+ :
S w(x, t) + w(x, t) + kw(x, t) + EIw (4) (x, t) = f0 (t)(x vt) Sg (5.1)
w(x, 0) = w0 ; w(x, 0) = v0
102 5. Etudes sur le modele 1D des voies
z
f0 (t)(x vt)
x +
wr
kt t semelles
wt mt traverses
kb b
ws ballast
ks s
sol
Ce modele simplifie considere la structure de couplage entre le rail, les semelles, les traverses
(blochets), le ballast et le sol comme un systeme 1D non-lineaire de poutres-masses-ressorts-
amortisseurs. Ce systeme est considere comme une structure qui est infinie et homogene en di-
rection horizontale. La figure 5.2 donne une illustration de ce modele ou on a utilise les hypotheses
suivantes :
Le rail est une poutre en flexion lineaire. Il peut alors etre modelise par une poutre dEuler-
Bernouilli. On designe par r , Er , Sr , Ir respectivement la masse volumique, le module dYoung,
laire de la section et le moment dinertie de flexion du rail.
Les semelles sont modelisees par un systeme uniforme de ressorts-amortisseurs k t , t .
Les traverses qui sont discretes dans la realite sont considerees comme un systeme lineaire continu
de masses reparties mt .
La couche de ballast ou la non-linearite sera introduite est modelisee par un systeme uniforme
et continu des barres verticales en negligeant les contraintes de cisaillement dans le ballast.
Les caracteristiques physiques de cette barre sont representees par la masse m b , la rigidite kb ,
lamortissement b . Dans le cas non-lineaire, la rigidite k b est supposee etre une fonction du type
unilateral de la deformation de la barre.
Le sol est suppose toujours lineaire et est modelise comme une fondation de Winkler dont la
rigidite et lamortissement sont respectivement k s et s .
La charge est une force ponctuelle qui se deplace en direction x avec une vitesse constante v.
Son amplitude est soit constante soit variable en fonction du temps.
Les inconnues a determiner dans le modele presente sont les deplacements, les vitesses et les
accelerations aux niveaux du rail, de la surface superieure du ballast (qui sont identiques a ceux
des traverses) et de la surface du sol. On note w r , wt et ws ces trois deplacements inconnus.
Afin dobtenir le systeme dequations dynamiques du probleme global, on va dabord ecrire
separement les equations dynamiques de chaque composante en fonction des forces de reaction
entre eux : le rail, les traverses et le ballast (Figure 5.3). En prenant en compte aussi la force
volumique, les equations dequilibre dynamiques secrivent comme suit :
5.1 Introduction 103
f0 (t)(x vt)
wr Fr
b , E b , b , Sb , L b
Fr
Ft
wt
wt
Ft Fs
ws
ws
Fs
Figure 5.3.
Pour le rail et les semelles : equation dynamique de la poutre dEuler-Bernouilli en flexion soumise
a une force f0 (t) et la force Fr distribuee par le systeme des ressorts :
r Sr (wr + g) + Er Ir wr(4) = Fr f0 (t)(x vt) (5.2)
Pour les traverses : equation dynamique des corps rigides avec les ressorts-amortisseurs :
mt (wt + g) + kt (wt wr ) + t (wt wr ) = Ft (5.3)
Pour le ballast : equation dynamique uniaxiale dune barre avec les parametres : b (masse
volumique), Eb (module delasticite), b (amortissement), Sb (aire de la section) et Lb (longueur).
En considerant cette barre comme un element barre de deux noeuds dont deux degres de liberte
sont wt et ws avec deux forces exterieures Ft , Fs , lequation dynamique en element fini secrit
(voir dans le chapitre 3 pour les matrices elementaires dune barre) :
mb 2 1 wt + g b 2 1 wt 1 1 wt Ft
+ + kb + =0
6 1 2 ws + g 6 1 2 ws 1 1 ws Fs
(5.4)
ou on a pose :
Eb Sb
mb = b Sb Lb ; b = b Lb ; kb = (5.5)
Lb
La notation g designe lacceleration de la gravite (g 9.81ms 2 ).
Il reste deux equations supplementaires qui permettent de calculer les forces dans les ressorts
des semelles kt (dont la deformation est (wr wt )) et du sol ks (dont la deformation est ws ) :
Fr = kt (wr wt ) t (wr wt ) (5.6)
Fs = ks ws s ws (5.7)
(x, t) R R+ :
mw(x, t) + w(x, t) + kw(x, t) + jw (4) (x, t) = f 0 (t)(x vt) g (5.9)
w(x, 0) = w 0 ; w(x, 0) = v 0
ou :
r Sr 0 0 t t 0
mb mb b b
m= 0 mt + = t t + (5.10)
3 6 3 6
mb mb b b
0 0 + s
6 3 6 3
kt kt 0 Er Ir 0 0
k = kt kt + kb kb j= 0 0 0 (5.11)
0 kb kb + k s 0 0 0
n mb mb o t
f = {f0 (t), 0, 0} tg = r Sr g, mt + g, g (5.12)
2 2
On trouve un modele a 3 inconnues du probleme general indique dans (4.4). Dans ce modele, on
suppose que la non-linearite napparat que dans la couche de ballast dont le comportement de
type unilateral est defini comme suit :
c
kb si wt ws
kb = t (5.13)
kb sinon
ou kbc , kbt sont deux rigidites differentes qui representent les resistances de la couche de ballast
en compression et en traction. De facon similaire, il faut definir deux valeurs differentes pour
lamortissement b en compression et en traction.
Remarque 5.1.1
Le vecteur w 0 de deplacement statique du a la force de gravite est constant pour tous les points x et est
evalue par :
w0 = k1 g (5.14)
=0 = 0.5
0.45 0.45
= 0.01
0.35 0.35
|w0/w0max|
0.25 = 0.25
0.25
= 0.25
0.2 0.2
= 0.5
= 0.75
0.15 = 0.75 0.15
= 0.5
0.1 =1 0.1 =1
=2 =2
0.05 0.05
0 0
0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4
=1 = 1.5
0.25
= 0.01
0.2 0.2
= 0.01
0.15 0.15
|w0/w0max|
|w0/w0max|
= 0.1
0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5
= 0.01 = 0.1
2 2 0.35
0.5 0.15
0.8 0.8
0.4
0.6 0.6 0.1
0.3
0.4 0.4
0.2
0.05
0.2 0.1 0.2
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
= 0.5 = 1.0
2 0.16 2 0.16
1.8 1.8
0.14 0.14
1.6 1.6
0.12 0.12
1.4 1.4
0.1 0.1
1.2 1.2
1 0.08 1 0.08
0
0.8 0.8
0.06 0.06
0.6 0.6
0.04 0.04
0.4 0.4
0.02 0.02
0.2 0.2
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
La valeur de la vitesse critique a un role tres important. On presente dans la figure (5.6)
les solutions des deplacements dus aux vitesses autour des valeurs critiques dans les 4 cas de
frequences precedents. Un coefficient damortissement = 0.1 est utilise. En passant la vitesse
5.3 Modele P1R 107
1 1
0
2 2
3 3
4 4
5 5 5
40 20 0 20 40 40 20 0 20 40 40 20 0 20 40
X (m) X (m) X (m)
1 0.5 0.5
deplacement vertical (m)
2 2 2
40 20 0 20 40 40 20 0 20 40 40 20 0 20 40
X (m) X (m) X (m)
1 1
0
0 0
1
1 1
2 2 2
40 20 0 20 40 40 20 0 20 40 40 20 0 20 40
X (m) X (m) X (m)
On va introduire la non-linearite dans cet exemple pour voir son influence sur une structure de
poutre due a la charge mobile. En supposant que le systeme de ressorts est plus faible en tension
pour representer le comportement de la couche de ballast, on va utiliser une rigidite de tension k t
qui est inferieure de celle de compression k c .
Il est evident que dans le cas non-lineaire, les valeurs des vitesses critiques sont changees. La
figure 5.7 presente deux courbes de deplacements qui nous permettent de determiner les vitesses
108 5. Etudes sur le modele 1D des voies
critiques dans deux cas : kt = 1%kc et kt = 50%kc . Seulement le cas dune charge constante a
ete traite. La vitesse critique diminue en fonction de la rigidite k t et vaut 80% de celle du cas
lineaire. Neanmoins, ce resultat na pas de sens general puisquil doit dependre surement aussi de
lamplitude de la charge.
= 0, k = 1%k = 0, k = 50%k
t c t c
1.8
1 = 0.1
= 0.1
1.6
1.4 = 0.25
0.8
1.2 = 0.25
|
|
0max
0max
1 0.6 = 0.5
|w /w
|w /w
0
0
= 0.5
0.8 = 1.0
= 1.0 0.4
0.6
0.4
0.2
0.2
0
0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4
v v
0 0
On presente dans la figure 5.8 deux exemples de cas ou la rigidite en tension des ressorts est
tres faible (kt = 1%kc ). Comme il ny a que la composante de force de gravite qui peut empecher
le mouvement de la poutre vers le haut, des que la deformation positive apparait dans les ressorts,
les deplacements deviennent tres importants.
3 = 0.5, = 0.5
4 = 0, = 2.0 x 10
x 10 3
20
lineaire lineaire
nonlineaire nonlineaire
15
2
deplacement vertical (m)
deplacement vertical (m)
10
1
0
0
5 1
200 150 100 50 0 50 100 150 200 50 0 50
X (m) X (m)
Dans la realite, la resistance en tension des ressorts nest pas vraiment nulle, grace aux forces
de cisaillement. La figure 5.9 montre les solutions non-lineaires comparees avec les calculs lineaires.
On impose une rigidite de tension qui vaut la moitie de celle de compression (k t = 50%kc ). La
structure est pre-contrainte sous la gravite avant le calcul dynamique, elle a un deplacement statique
w0 < 0. La condition pour verifier si le ressort en un point est en tension est realisee simplement
par w(x) + w0 >? 0. Les valeurs de vitesse (v) et de frequence () utilisees dans ces calculs sont
5.3 Modele P1R 109
4 = 0, = 0.5 4 = 0, = 2.0
x 10 x 10
0 4
lineaire
nonlineaire
3
2 0
3 2
3
lineaire
nonlineaire
4 4
10 5 0 5 10 40 30 20 10 0 10 20 30 40
X (m) X (m)
4 = 0.5, = 0.5 4 = 0, = 1.2
x 10 x 10
4 1
lineaire lineaire
nonlineaire nonlineaire
0.5
2
0
0
0.5
2
1
4
1.5
6
2
8 2.5
40 30 20 10 0 10 20 30 40 40 30 20 10 0 10 20 30 40
X (m) X (m)
2 1
1
0
1
1
2 2
40 30 20 10 0 10 20 30 40 40 30 20 10 0 10 20 30 40
X (m) X (m)
2
1
0 1
1
0
2
1
3
4 2
40 30 20 10 0 10 20 30 40 40 30 20 10 0 10 20 30 40
X (m) X (m)
Figure 5.9. Influence de la vitesse dune charge harmonique pour le rail pose
sur une fondation non-lineaire de Winkler (k t = 50% kc )
110 5. Etudes sur le modele 1D des voies
toujours calculees avec les valeurs critiques qui correspondent a la rigidite en compression (k c ).
On constate que la presence de la non-linearite ne change presque pas le mecanisme de pro-
pagation dondes. Par contre, elle augmente lamplitude du deplacement. La longueur donde est
aussi plus grande dans une structure non-lineaire.
On a vu dans la section 2.5 une formule analytique de force exercee par le rail sur un blochet
qui est proposee par Alaoui et Naciri [4] en utilisant des calculs statiques. On va chercher dans
cette partie les surcharges dynamiques dues a la vitesse.
On va determiner dabord la rigidite equivalente du sol k qui donne la meme force que celle
donnee dans [4] dans le cas statique. Rappelons que le signal de force 2.89 du a deux forces de Q/2
a la distance L = 3m secrit sous la forme suivante :
QY h ( V ta )2 V taL 2
i
F3d (t, Q, V ) = X d + X( d ) (5.15)
2
et si lon remplace le sol par un systeme de ressorts, cette force peut etre aussi determinee analy-
tiquement (voir [54]) en faisant attention que la distance entre les traverses est egale a d :
Qd |V t a| |V t a| |V t a|
F1d (t, Q, V ) = exp cos + sin +
4
|V t a L| |V t a L| |V t a L|
+ exp cos + sin
(5.16)
q
ou = 4 4EI k . On peut en deduire directement la condition pour que les maximums des deux
fonctions F3d et F1d soient les memes :
QY Qd
= (5.17)
2 4
qui donne la rigidite equivalente k :
4
2Y
k = 4EI (5.18)
d
On presente dans la figure 5.10 la comparaison des forces sur un blochet en fonction du temps
dans les deux modeles 3D et 1D. Deux rigidites du sol : E s = 100 M P a et Es = 30 M P a sont
utilisees qui donnent respectivement deux valeurs de k : k = 1.0310 8 N/m et k = 3.17107 N/m.
On constate que la fondation de Winkler represente assez bien la fondation 3D si le sol est rigide.
La figure 5.11 presente la comparaison entre les charges statiques et dynamiques. Les calculs
sont effectues avec une vitesse de 200 m.s 1 en tenant compte de la masse des traverses. Les
resultats montrent des differences sensibles dues aux surcharges dynamiques, surtout quand le sol
est souple (le cas Es = 30 M P a).
Le modele 1D complexe (rail, semelles, blochets, ballast, sol) sera etudie dans cette partie. Les
caracteristiques de la structure sont donnees dans le tableau 5.1. Ces donnees, en fait, representent
la structure totale de voie comme un systeme de poutre-ressort-amortisseur comme suit :
5.4 Modele P3R 111
2 1.5
force (N)
force (N)
1.5 1
1
0.5
0.5
0
0
0.5 0.5
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
temps (s) temps (s)
force (N)
2 1.5
1.5
1
1
0.5
0.5
0
0
0.5 0.5
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
temps (s) temps (s)
La couche de ballast dans la realite est tres absorbante. On prend ici pour elle une valeur damor-
tissement b tres importante (b = 1).
On regarde dans cette partie les reponses de ce systeme en fonction des types de charges
(vitesse, frequence) et aussi du comportement de la structure (raideurs, amortissements ...) dans
les cas lineaires et non-lineaires.
Dans le calcul statique des voies ferrees, chaque essieu supporte une charge P = 170 kN . Pour
les exemples ici, on va prendre :
Cas statique : f = P .
P
Cas dynamique : f (t) = P + 2 sin t.
Influence des vitesses et des frequences. Le premier exemple est realise en supposant que
les rigidites du ballast et du sol sont k b = 16 108 kN.m1 (pour un module dYoung du ballast
de 200 M P a) et ks = 5 107 kN.m1 . Dans tous les calculs, les deplacements statiques dus aux
poids de la structure sont pris en compte. Pour cette structure, ils valent : w r = 0.366 mm,
wt = 0.3625 mm et ws = 0.355 mm.
La structure est calculee avec 2 vitesses differentes v = 100, 200 m.s 1 et 3 frequences differentes
f0 = 0, 5, 10 Hz. Pour les cas dynamiques, les calculs sont effectues pendant 4 periodes afin darriver
2
au regime stationnaire (des pas de temps qui valent dt = 0.05 sont utilises pour trouver des
resultats stables dans le schema de Newmark).
Sur le deplacement. Les resultats des deplacements aux trois niveaux w r , wt et ws par elements
finis sont presentes dans la figure 5.12. Les solutions sont tracees aux derniers pas de temps (quand
f (t) = P ).
On constate que deux valeurs de deplacement aux niveaux des traverses w t et du sol ws sont tres
proches car la couche de ballast est tres rigide par rapport aux semelles et au sol. Les deplacements
qui decroissent tres lentement derriere la force montrent leffet de lamortissement important dans
la structure. Ainsi, le deplacement maximum ne change pas beaucoup en fonction de la vitesse et
de la frequence grace aux amortissements.
Quand la charge devient harmonique, la vibration de la partie X < 0 est influencee nettement
par la frequence f0 . Au contraire, la partie X > 0 ne change presque pas.
Sur lacceleration et la deformation de la couche de ballast. Pour la couche de ballast, il est aussi
interessant de considerer les deux cas dangereux suivants :
- quand lacceleration verticale positive (vers en haut) dans la couche de ballast depasse la gravite :
les cailloux peuvent sauter. On verifie la condition : t = wt > g pour cette situation.
- quand la couche de ballast nest plus en compression, une discontinuite va apparatre dans cette
couche. On verifie la condition wt ws > 0 pour cette situation.
5.4 Modele P3R 113
4
v = 100 ms1, f = 0 Hz 4
v = 200 ms1, f = 0 Hz
0 0
x 10 x 10
2 3
3
4
4
deplacement vertical (m)
6
6
7
7
8
9 8
10
rail 9 rail
11 traverse traverse
sol sol
10
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
1 1
4
v = 100 ms , f0 = 10 Hz 4
v = 200 ms , f0 = 10 Hz
x 10 x 10
2 3
3
4
4
deplacement vertical (m)
6
6
7
7
8
9 8
10
rail 9 rail
11 traverse traverse
sol sol
10
20 15 10 5 0 5 10 15 20 40 30 20 10 0 10 20 30 40
X (m) X (m)
fo = 0 Hz fo = 10 Hz
10 20
1 1
v = 50 ms v = 50 ms
8 v = 100 ms
1
v = 100 ms
1
15
6 v = 200 ms1 v = 200 ms1
4 10
2
5
(m.s )
(m.s )
2
0
0
t
4 5
6
10
8
10 15
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
4 fo = 0 Hz 4 fo = 10 Hz
x 10 x 10
0 0
1
deformation du ballast
deformation du ballast
1 2
3
v = 50 ms1 v = 50 ms1
1 1
v = 100 ms v = 100 ms
v = 200 ms1 v = 200 ms1
2 4
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
On evalue a chaque pas de temps deux fonctions t et (wt ws )/hb (hb = 30 cm est la hauteur
de la couche de ballast) pour chercher linstant ou elles ont des valeurs maximales. Dans la figure
5.13, les courbes de ces 2 fonctions a cet instant sont tracees. Dans le cas dune charge constante, la
valeur maximale de vitesse de la charge pour le critere w t < g vaut environ 200 ms1 . En revanche,
cette valeur devient inferieure a 100 m.s 1 quand la charge est harmonique.
Les courbes de (wt ws )/hb dans la figure 5.13 montrent que la couche de ballast reste toujours
en compression. La non-linearite napparat pas dans ce probleme.
Effet de la rigidite du sol. On prend dans cet exemple un sol qui est plus souple et dont
la rigidite ks = 2 107 N m1 . La rigidite du ballast est toujours k b = 1.6 109 N m1 . Les
deplacements sont evidemment plus importants. Par contre, les accelerations (figure 5.17) sont
moins importantes que celles dans le cas ou le sol est plus rigide (figure 5.13).
Effet de deux charges. On considere dans ce calcul le probleme de deux charges qui representent
deux essieux dont la distance entre eux est de 3 m. On utilise donc deux forces identiques qui
sappliquent aux points X = 1.5 m et X = 1.5 m. Les memes calculs sont realises (avec le ballast
rigide).
Les deplacements verticaux maximums (figure 5.18) dus a deux forces ne changent pas vis a
vis de ceux dus a une seule force. La meme remarque sapplique pour lacceleration.
La figure 5.20 presente la force au niveau du haut de la couche de ballast pour avoir une idee
de la surcharge dynamique a cause de la vitesse et de la frequence.
Remarque 5.4.1
En regardant les exemples presentes, on peut faire les remarques suivantes sur les resultats :
Le deplacement est plus important avec une structure plus souple.
Les accelerations verticales en haut de la couche ballast sont plus importantes si le ballast est plus souple.
En revanche, elles ne sont pas beaucoup imfluencees par un sol plus souple.
5.4 Modele P3R 115
4
v = 100 ms1, f = 0 Hz 4
v = 200 ms1, f = 0 Hz
0 0
x 10 x 10
2 3
4
4
5
deplacement vertical (m)
8 7
10 8
9
12
10
14 rail rail
traverse 11 traverse
sol sol
16
20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20
X (m) X (m)
4
v = 100 ms1, f = 10 Hz 4
v = 200 ms1, f = 10 Hz
0 0
x 10 x 10
2 3
4
4
5
deplacement vertical (m)
7
8
8
10 9
10
12 rail rail
traverse 11 traverse
sol sol
14
20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20
X (m) X (m)
fo = 0 Hz fo = 10 Hz
25 30
1 1
v = 50 ms v = 50 ms
1 1
20 v = 100 ms v = 100 ms
v = 200 ms1 20 v = 200 ms1
15
10 10
(m.s )
(m.s )
5
2
0
0
t
5 10
10
20
15
20 30
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
4 fo = 0 Hz 3 fo = 10 Hz
x 10 x 10
0 0
2
deformation du ballast
deformation du ballast
0.5
3
5
1
6
v = 50 ms1 v = 50 ms1
7 v = 100 ms
1
v = 100 ms1
v = 200 ms1 v = 200 ms1
8 1.5
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
1
v = 100 ms , f = 0 Hz 1
x 10
3 0
3 v = 200 ms , f = 0 Hz
0
x 10
0.8
0.9
1
deplacement vertical (m)
1.1
1.2
1.5
1.3
rail
1.4 rail
traverse
traverse
sol
sol
20 15 10 5 0 5 10 15 20 1.5
X (m) 20 15 10 5 0 5 10 15 20
X (m)
1 1
3 v = 100 ms , f = 10 Hz 3 v = 200 ms , f = 10 Hz
0 0
x 10 x 10
0.8 0.5
1
deplacement vertical (m)
1.4
1.6 1.5
fo = 0 Hz fo = 10 Hz
10 20
1 1
v = 50 ms v = 50 ms
8 v = 100 ms
1
v = 100 ms
1
15
6 v = 200 ms1 v = 200 ms1
4 10
2
5
(m.s )
(m.s )
2
0
0
t
4 5
6
10
8
10 15
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
4 fo = 0 Hz 4 fo = 10 Hz
x 10 x 10
0 0
0.5
deformation du ballast
deformation du ballast
1 1.5
1 1
4 v = 100 ms , f = 0 Hz 4 v = 200 ms , f = 0 Hz
0 0
x 10 x 10
2 2
4 4
deplacement vertical (m)
8
8
10
10
12
12
14 rail rail
traverse 14 traverse
sol sol
16
20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20
X (m) X (m)
1
1 v = 200 ms , f = 10 Hz
4
v = 100 ms , f = 10 Hz x 10
4 0
0
x 10
2
10
10
12 rail
rail
traverse
traverse
sol
sol
14 20 15 10 5 0 5 10 15 20
20 15 10 5 0 5 10 15 20 X (m)
X (m)
fo = 0 Hz fo = 10 Hz
15 12
v = 50 ms1 v = 50 ms
1
5 4
t (m.s2)
(m.s2)
2
t
0 0
5 4
10 8
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
4 fo = 0 Hz 4 fo = 10 Hz
x 10 x 10
0 0
0.5
1
deformation du ballast
deformation du ballast
1.5
3
2
v = 50 ms1 v = 50 ms1
v = 100 ms1 v = 100 ms1
v = 200 ms1 v = 200 ms1
2.5 4
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
Figure 5.19. Acceleration t = wt et deformation du ballast (wt ws )/hb : cas de deux essieux
118 5. Etudes sur le modele 1D des voies
0
0
2
4
5
force (N)
force (N)
6
10
8
10
15
12
statique statique
dynamique dynamique
14 20
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
5.5.1 Modele
On va considerer dans cette partie le cas ou le systeme de ressorts nest plus uniforme mais
periodique pour representer les traverses des voies ferrees. On va utiliser un modele de calcul
presente dans la figure 5.21 ou la poutre est supportee par des morceaux des ressorts uniformes
de largeur 30 cm qui sont mis chaque 60 cm 1 . En consequence, la rigidite et lamortissement des
ressorts sont pris egaux a 2k et 2, respectivement.
z f0 (t)(x vt)
x
+
2k 30 cm 30 cm
2
Un calcul transitoire dans le repere fixe doit etre envisage pour traiter ce probleme. Le probleme
peut etre resolu par une meme procedure que celle presentee dans le cas uniaxial (section 3.3).
On cherche sil y a des differences de reponse dynamique dans la structure avec les modeles
continus et discrets. Dans les calculs en dessous, on utilise un maillage de 800 elements identiques
dans un domaine dont la taille de chaque element vaut h = 7.5 cm. Le pas de temps de calcul est
pris egal a t = h/v ou v est la vitesse de la charge.
Deux cas de rigidite de la fondation sont consideres. Le premier calcul est avec les donnees
du tableau 4.1 dont la rigidite de la fondation k vaut 1.6 10 7 N/m. Dans le deuxieme cas, on
1
Dans la realite, les largeurs des traverses sont de 29 cm et la distance entre elles est de 60 cm
5.5 Probleme dappuis discrets 119
x 10
4 = 0, = 0.5 4 = 0, = 0.5
1 x 10
3
2
0
1
2
2
3
3
uniforme uniforme
periodique periodique
4 4
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
x 10
4 = 0.5, = 0.5 4 = 0.5, = 1.5
2 x 10
4
1
deplacement vertical (m)
0
3
4
uniforme uniforme
periodique periodique
5 2
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
0
deplacement vertical (m)
deplacement vertical (m)
2
0
4
10
1
12
uniforme uniforme
14 periodique periodique
5 0 5 5 0 5
X (m) X (m)
x 10
5 = 0.5, = 0.5 x 10
6 = 0.5, = 1.5
1 8
4
0.5
deplacement vertical (m)
deplacement vertical (m)
0
0
2
0.5
6
8
uniforme uniforme
periodique periodique
1 10
5 0 5 5 0 5
X (m) X (m)
calcule avec une rigidite k 100 fois plus grande, i.e. k = 1.6 10 9 N/m. La structure est calculee
en subsonique et en supersonique, avec une amplitude de charge soit constante soit harmonique.
On peut constater en regardant la figure 5.22 que, dans le cas ou la fondation est souple, la
solution obtenue est la meme avec un modele dappuis continus ou avec un modele dappuis discrets.
En revanche, les differences deviennent visibles quand la rigidite des appuis est importante. En effet,
quand la fondation est souple, la deflexion locale de la poutre entre deux appuis est tres petite par
rapport a la deflexion globale. Ce rapport augmente en fonction de la rigidite des appuis.
5.6 Conclusion
Les etudes sur le modele unidimensionnel des voies ferrees sont effectuees avec deux modeles :
P1R et P3R. A partir de ces resultats, on peut deduire les conclusions suivantes :
La vitesse critique a un role tres important. Des que la vitesse de la charge depasse la vitesse
critique, la nature de la solution change completement. Dans le cas non-lineaire, la valeur de la
vitesse critique depend aussi de lamplitude de la charge.
Les deplacements sont moins importants dans le cas ou le ballast et le sol sont plus rigides.
La rigidite faible de la couche de ballast peut creer des accelerations verticales importantes.
On na pas trouve de probleme ou il existe des ondes de choc comme dans le cas uniaxial (on
note aussi que la non-linearite cette fois-ci est introduite dans les ressorts, pas dans la structure
comme pour le cas uniaxial).
Comme la distance entre les blochets nest pas grande par rapport a leur taille, leffet de
periodicite des appuis discrets nest pas tres important dans les calculs realises.
Troisieme partie
Modelisation du comportement 3D
du ballast
6.1 Introduction
6.2.1 Motivation
On sinteresse a une loi de comportement macroscopique des materiaux granulaires qui adapte
leffet non-tension pour ce type de materiau. Considerons dabord le cas unidimensionnel ou il y
a un deplacement impose (L) sur une chane de grains de longueur L (Figure 6.1). On definit la
deformation par = L L et la contrainte par la valeur moyenne des forces de contact entre les
grains. Cette chane de grains verifie une loi de comportement du type unilaterale suivante :
= k si <0
= 0 si 0
ou k est la rigidite moyenne de contact des grains. La figure 6.1 presente une courbe de compor-
tement unilateral dans le cas ou k est lineaire. Cette courbe montre le caractere de non rigidite en
tension des materiaux granulaires.
On integre cette propriete dans le cas plus general du milieu granulaire en 3D. A partir du fait
quil ny pas de resistance entre deux grains en traction, autrement dit, quon ne peut pas appliquer
une contrainte positive sur les materiaux granulaires, on propose un modele de comportement
macroscopique que lon appelle modele non-tension. Dans ce modele, la loi de comportement doit
etre introduite telle que la contrainte normale sur nimporte quelle section dans le domaine considere
ne soit jamais positive ( ij ni nj 0 ou et n sont le tenseur des contraintes et le vecteur unite
normal de la section consideree).
124 6. Modelisation du comportement 3D du ballast
Dans la litterature, le modele non-tension a ete propose depuis quelques annees, surtout dans
les travaux sur la modelisation des structures masonry-like. Pazenca, Polizzotto [91] et Del Piero
[29] ont discute de la loi de comportement et aussi de la compatibilite de la charge exterieure pour
ce modele ou la deformation est decomposee en deux parties : lune elastique et lautre anelastique.
Les hypotheses fondamentales de non-resistence en tension et dexistence de la densite denergie
elastique demandent que le tenseur des contraintes soit defini negatif et que la contrainte et la
deformation anelastique soient orthogonales. Lhypothese de petite deformation est toujours uti-
lisee. Suite a cette base theorique, les aspects numeriques sont aussi etudies dans le developpement
dune formulation variationnelle et dans la mise en oeuvre dans la methode des elements finis
[5][26][77].
La loi de comportement du modele non-tension des materiaux granulaires proposee ici est
derivee directement de la densite denergie de deformation qui sexprime en fonction des deformations
principales. Cette fonction sannule en fonction du signe des deformations principales et aussi de
la dilatation volumique. Le probleme est traite comme dans le cas des materiaux hyperelastiques.
La loi de comportement secrit dans la base principale des tenseurs des deformations a laide de la
decomposition spectrale et ensuite est integree dans la procedure de la methode des elements finis
non-lineaires.
Dans les calculs suivants, on appelle S 3 lespace des tenseurs du 2e ordre de dimension 3. Les
notations S3+ , S3 indiquent lespace des tenseurs definis positifs ou negatifs dans S 3 .
A. Hypotheses
(i) Soit V un volume unite dans une structure de materiaux granulaires. En supposant que le
nombre de grains dans V est suffisamment grand tel que les torsions locales peuvent etre
negligees, la direction des forces de contact se distribue de facon isotrope. Le domaine peut
etre considere comme un milieu isotrope et homogene.
(ii) La deformation est infinitesimale.
(iii) Le tenseur des contraintes est toujours defini negatif.
Pour decrire le comportement de materiaux elastiques normaux, il suffit dintroduire une fonc-
tion denergie des deformations qui est une forme quadratique positive convexe de la deformation
W () : S3 R+ et sexprime par :
ou C est un tenseur du 4e ordre, dont les composantes sont les coefficients elastiques et , S 3
sont les tenseurs des contraintes et des deformations. Par definition, le tenseur des contraintes est
obtenu a partir de la derivee de lenergie W :
W
= (6.2)
Comme doit etre toujours defini negatif, il est preferable decrire W en fonction des valeurs
principales de deformation. Pour sadapter aux materiaux granulaires, on suppose quil existe une
energie non nulle dans une direction n si soit la deformation dans cette direction ( ij ni nj ) soit la
dilatation volumique est negative. On propose alors une fonction W qui a la meme forme que dans
le cas elastique classique :
3
1 X
W () = [(1 + 2 + 3 ) ]2 + [
a]
2
(6.3)
2 a=1
1
= ( ||) (6.4)
2
Dans le cadre de ce chapitre, on considere deux cas differents qui dependent des hypotheses
utilisees pour et :
Elasticite lineaire. et sont deux coefficients constants (deux coefficients elastiques de Lame).
Elasticite non-lineaire. est une fonction de tr : = (tr) et est une fonction qui est definie
telle que : [ 2 2
a ] = (a ) [a ] .
On va dans la suite exprimer les derivees de W dans ces deux cas lineaire et non-lineaire. Ces
derivees vont servir aux calculs suivants.
= H() ; = (6.5)
ou H() est la fonction de Heaviside qui est definie par :
1 si x0
H(x) = (6.6)
0 si x<0
W
= (tr) + 2
a
a
(6.7)
2 W
= H(tr) + 2 H(a )ab
a b
W
= 1 (tr) + 21a
a
a
(6.8)
2 W
= 2 H(tr) + 22a H(a )ab
a b
avec :
1
1 = + 0 tr (6.9)
2
1
1a = a + 0a a (6.10)
2
1
2 = + 20 (tr) + 00 (tr)2 (6.11)
2
0 1 00 2
2a = a + 2a a + a a (6.12)
2
ou on note : a = (a ).
Demonstration.
En utilisant (6.5) et notant que = , la derivee de la fonction de densite denergie W (6.3)
peut etre calculee comme suit :
W (tr) (tr) 1 0 2 2
= (tr) + (tr) + 2a a
a
+ 0a
a
a (tr) a 2 a
1 1
= [ + 0 tr ]tr + 2 [a + 0a a ] a
2 2
1 1
= [ + 0 tr]tr + 2 [a + 0a a ]
a (6.13)
2 2
La derivee du deuxieme ordre est :
W 1 0 tr tr 1 1
= [ + tr] + [0 + 0 + 00 tr] tr
b a 2 tr b 2 2
1 0 a 1 1
+ 2 [a + a a ] a + 2 [0a + 0a + 00a a ]
a (6.14)
2 a b 2 2
2 W 1 1
= [ + 20 (tr) + 00 (tr)2 ]H(tr) + 2 [a + 20a a + 00a 2a ] H(a )ab (6.15)
a b 2 2
Remarque 6.2.1
1. La fonction de densite denergie proposee pour le cas non-lineaire assure toujours lisotropie du comporte-
ment du materiau. Cela peut etre demontre en verifiant la condition pour quune fonction soit isotropique
qui demande quon ait Q orthogonal direct 1 :
1
Un tenseur ' est appele orthogonal direct sil verifie : '(' t
= ' '
t
'
= 1 et det( ) = 1
6.3 Methode de resolution numerique 127
2. Dans le domaine des petites deformations, on peut ecrire les deux coefficients delasticite et i comme
des polynomes en fonction de tr et i . En fait, les donnees experimentales nous donnent toujours la
relation de contrainte-deformation, cest a dire les fonctions 1 et 1a . Supposons que 1 soit un polynome
du 4e ordre en tr :
6.3.1 Rappel de la methode des elements finis pour le probleme elastique non-
lineaire
La loi de comportement des materiaux granulaires a ete decrite comme un materiau non-lineaire
elastique. On va presenter la procedure numerique generale de la methode des elements finis dans
laquelle on peut integrer la loi de comportement proposee dans les calculs de structures comportant
des materiaux granulaires.
Soit R3 un domaine continu dont la frontiere = u f . Les conditions aux limites sont
les forces appliquees sur f et le deplacement nul impose sur u qui secrivent :
Appelons u(x) le vecteur deplacement en un point x . Le Principe des Travaux Virtuels (PTV)
sur une perturbation virtuelle u statique secrit :
Z Z
G(u, u) =
(x)T d
f(x)
t u d = 0
f
ad
x ; u C (6.22)
ou C ad est lespace admissible du deplacement virtuel. Dans cette equation, on a ecrit les tenseurs
des contraintes et des deformations sous la forme de vecteurs pour simplifier les calculs dans
la suite.
En utilisant la discretisation de la methode des elements finis, le travail virtuel peut se reecrire
en fonction du vecteur des deplacements aux noeuds U :
ou R(U ) est appele le vecteur residu. Le PTV impose R(U ) = 0 et il faut resoudre cette equation
pour trouver la solution U . Les algorithmes iteratifs sont utilises afin de trouver la solution U
dans le cas non-lineaire. On utilise ici la methode de Newton - Raphson qui est basee sur le
developpement de Taylor au 1e ordre de R(U ) a partir dune solution actuelle U (k) .
(k+1) (k)
R
) )
R R + U k+1 = 0 (6.24)
U = (k)
128 6. Modelisation du comportement 3D du ballast
ou le residu R(U )(k) est determine en fonction des deplacements a lincrement precedent U (k) :
Z
(k)
R = B t (k) d F (k) (6.26)
on rappelle que la fonction dinterpolation nous permet de calculer la deformation a partir des
deplacements aux noeuds = BU . En supposant que la force exterieure ne depend pas du
(k)
deplacement, la matrice de rigidite tangente K T est determinee par :
Z (k) Z
(k) R(k)
KT = = B t
d = B t D (k) B d (6.27)
U U
*
et D (k) = (k) definit le module tangent du materiau. D (k) est la presentation matricielle du
tenseur tangent du 4e ordre qui est defini par 2 W . *
Les deux termes a la droite de lequation (6.25) sont determines explicitement en fonction de
(k)
U (k) . On peut alors calculer lincrement U (k+1) a partir de linverse de K T et puis calculer le
vecteur des deplacements a literation (k + 1) :
La solution finale du deplacement est calculee de facon iterative jusqua la convergence qui est
verifiee par la condition :
|U (k+1) |
< (6.29)
|U (k+1) |
ou est la tolerance.
On a presente ci-dessus la methode generale des elements finis pour resoudre le probleme
elastique non-lineaire. Le probleme restant a resoudre est comment calculer le residu et le module
tangent avec le comportement non-tension a chaque iteration (k) (cest a dire comment calculer
et C). La partie suivante presente une methode de decomposition spectrale qui permet dexpliciter
et C a partir dun tenseur des deformations donne.
3
X
= a n(a) n(a) kn(a) k = 1 (6.30)
a=1
ou n(a) (a = 1, 2, 3) sont les vecteurs principaux de et a sont les valeurs principales associees 2 .
Ces valeurs principales sont en fait les racines du polynome caracteristique :
ou les coefficients I1 , I2 , I3 sont appeles les invariants du tenseur et sont definis par :
1
I1 = : 1 ; I2 = [( : 1)2 2 : 1] ; I3 = det [] (6.32)
2
ou 1 est le tenseur identite du 2e ordre. Si on appelle 1 , 2 , 3 les trois valeurs principales de ,
I1 , I2 , I3 sexpriment par :
I1 = 1 + 2 + 3 (6.33)
I2 = 1 2 + 2 3 + 3 1 (6.34)
I3 = 1 2 3 (6.35)
Les vecteurs {n(a) , a = 1, 2, 3} forment une base orthonormee et verifient les quelques proprietes
utiles suivantes :
3
X
(i) n(a) n(a) = 1 (6.36)
a=1
(a)
(ii) n .n(b) = ab (6.37)
Remarque 6.3.1
Soient 1 , 2 , 3 les valeurs principales du tenseur . Notons {a, b, c} une permutation de {1, 2, 3}, on distingue
les 3 cas suivants :
(i) a 6= b 6= c : Les vecteurs principaux sont independants et la base orthogonale {n(a) , n(b) , n(c) } est
specifiee dans R3 en fonction des valeurs de a , b , c (figure 6.2a).
(ii) a = b 6= c : n(a) et n(b) sont dependants et en utilisant (6.36) dans (6.30), le tenseur peut secrire
dans ce cas :
= a 1 n(c) n(c) + c n(c) n(c) (6.38)
Le seul vecteur qui est specifie est n(c) . Les deux restants peuvent tourner librement dans le plan
perpendiculaire au vecteur n(c) (figure 6.2b).
(iii) a = b = c : on peut choisir nimporte quelle base orthogonale dans R3 pour {n(a) , n(b) , n(c) } (figure
6.2c) et le tenseur sexprime par :
= a 1 (6.39)
b
c c b c
a
a
a
A
, Proposition
P 6.1
Soit = 3a=1 a n(a) n(a) la decomposition spectrale de , si a est une valeur principale unique
de on peut ecrire lexpression explicite du tenseur du 2 e ordre m(a) := n(a) n(a) en fonction de
130 6. Modelisation du comportement 3D du ballast
de la facon suivante :
avec Da = 2a I1 + I3 2
a
Demonstration.
P3 P3
Par definition, on a a=1 m(a) = a=1 n(a) n(a) = 1. Multiplions m(a) par le produit (a b )(a c ) :
Remarque 6.3.2
La proposition (6.40) est validee si et seulement si Da 6= 0. En notant que Da peut autrement secrire
par Da = (a b )(a c )/a , la condition dexistence de m(a) est verifiee si a est une racine simple de
lequation p() = 0.
, Proposition 6.2
Soit a une valeur principale de . A partir de (6.31), a est une fonction de , a = a () et ses
derivees sexpriment par :
a
= n(a) n(a) m(a) (6.42)
Demonstration.
La derivee de la formule spectrale de nous donne :
3
X
d = [db n(b) n(b) + b dn(b) n(b) + b n(b) dn(b) ] (6.43)
b=1
en contractant cette equation avec n(a) n(a) et en utilisant le fait que kn(a) k = 1 n.dn = 0,
on obtient :
Remarque 6.3.3
Suite a cette proposition, on peut aussi determiner les derivees des invariants I 1 , I2 , I3 en fonction de dans
le cas 1 6= 2 6= 3 comme suit :
I1 I2 I3
=1 ; = I1 1 ; = I3 1 (6.45)
6.3 Methode de resolution numerique 131
, Proposition 6.3
La derivee de m(a) en fonction de sexprime explicitement par :
*
m(a) =
1 h
Da
I 1 1 + I3 1
a
1 0
1 Da m(a) m(a)
*
i
(a) (a)
+ I3 1
a I 1 + m 1 + 1 m (6.46)
I3 2
a
m(a) 1 + 1 m(a)
Da
ou Da = 2a I1 + I3 2
a ; Da0 = 2 2I3 3
a ; I est le tenseur identite du quatrieme ordre :
1
Iijkl = (ik jl + il jk ) (6.47)
2
ou ij est le delta de Kronecker ; la notation I * 1 designe un tenseur du quatrieme ordre qui est
defini par :
(I * 1
1
)ijkl = (1 1 + 1
2 ik jl
1
il jk ) (6.48)
Demonstration.
A partir de (6.40), on a : Da m(a) = (I1 a ) 1 + I3 1
a
1
et pour montrer (6.46), il suffit dexprimer
(a)
[Da m ] soit :
autrement, on a :
(a) Si a 6= b 6= c
Par definition, la contrainte est la derivee au 1 e ordre de la fonction de densite denergie (6.2).
Comme a 6= b 6= c , on peut decomposer cette derivee dans la base principale :
3
W X W a
= = (6.51)
a
a=1
*
sachant que a = m(a) (6.42), on obtient :
3
X
= W ,a n(a) n(a) (6.52)
a=1
P
Cela represente la decomposition spectrale de : = 3a=1 a n(a) n(a) avec a = W ,a qui sont
les valeurs principales de . On voit bien que les directions principales du tenseur des contraintes
sont identiques a celles du tenseur de deformation.
On derive encore une fois la fonction energie de deformation afin de calculer le module tangent :
3
X
C = ( W ,a m(a) )
a=1
3 3
X W ,a X m(a)
= m(a) + W ,a (6.53)
a=1 a=1
en developpant cette expression comme dans (6.51), on peut exprimer le module tangent :
*
X 3
3 X 3
X
C= W ,ab m(a) m(b) + W ,a m(a) (6.54)
a=1 b=1 a=1
(b) Si a = b 6= c
Si le tenseur a deux (et seulement deux) valeurs principales identiques, il ny a que le vecteur
principal dans lautre direction qui peut etre determine. Comme on a vu dans (6.38), la contrainte
peut aussi sexprimer en fonction de la troisieme direction principale :
a b
= = W ,aa (6.60)
a b
a b
= = W ,ab (6.61)
b a
a b c c
= = = = W ,ac (6.62)
c c a b
c
= W ,cc (6.63)
c
3 3 3
X X X
I = = i m(i) = m(i) m(i) + i m(i)
i=1 i=1 i=1
X
= m(i) m(i) + m(c) m(c) + a (m(a) + m(b) ) + c m(c)
i=a,b
X
= m(i) m(i) + m(c) m(c) + a (1 m(c) ) + c m(c)
i=a,b
X
= m(i) m(i) + m(c) m(c) + (c a ) m(c) (6.66)
i=a,b
qui donne :
X
m(i) m(i) = I m(c) m(c) (c a ) m(c) (6.67)
i=a,b
on en deduit :
(c) Si a = b = c
= W ,a 1
(6.73)
C = (W ,aa W ,ab )I + W ,ab 1 1
En resume, pour un materiau dont la densite denergie est une fonction des deformations prin-
cipales, a chaque etat de deformation , on peut calculer explicitement et C. La procedure de
calcul est presentee dans le tableau 6.1.
Remarque 6.3.4
Dans les cas ou la difference entre 2 deformations a , b est petite, cela induit des singularites numeriques
dans le calcul de m(a) , m(b) si on utilise toujours la formule (6.40) . Dans la mise en uvre numerique, on
va considerer a , b comme des valeurs identiques si (|a |, |b |) verifient la condition suivante :
|a b | TOL max (6.74)
3
avec max = max(|a |, |b |, |c |) et TOL = 10
6.4 Identification du comportement statique du ballast 135
On va presenter dans cette partie des travaux experimentaux afin de trouver le comportement
statique du micro ballast. Les essais uniaxiaux sont realises avec du ballast mis dans des tubes
cylindriques. Les charges sont appliquees en plusieurs cycles afin datteindre un etat stable pour le
ballast.
capteur de dplacement
jauge transversale
jauge longitudinale
capteur de
dplacement
D1
tube de PVC
J2 J1
BALLAST
D2
Les caracteristiques mecaniques du PVC utilise sont determinees simplement par un essai de
compression du tube seul. Il nous donne une courbe lineaire qui est la relation entre la deformation
axiale et la force verticale appliquee. La pente de cette courbe nous permet de determiner la rigidite
verticale du tube EverP V C = 1.23 GP a (E P V C = / ). Le coefficient de Poisson du PVC est
ver 33 33
obtenu en mesurant la deformation laterale et vaut = 0.4.
Le ballast utilise dans cet essai est a lechelle 1/3 et resulte dun melange de cailloux avec les
3 differentes coupures suivantes :
20 % pour la classe granulometrique 6-10 mm
30 % pour la classe granulometrique 10-14 mm
50 % pour la classe granulometrique 14-20 mm
Comme le comportement du ballast depend beaucoup de sa mise en uvre, les cailloux mis
dans le tube de PVC sont dabord pre-tasses en utilisant une machine de tassement.
La deformation laterale du ballast est recuperee a partir de la deformation du tube de PVC.
Quand on applique une force verticale sur le ballast, elle va causer deux deformations dans le tube
de PVC : la deformation radiale due a la pression interne et la deformation verticale due a la force de
frottement entre les cailloux et le tube. Pour mesurer ces deformations, deux jauges longitudinales
(J1) et deux jauges transversales (J2) sont donc collees symetriquement sur la surface exterieure
du tube PVC au niveau du milieu de la hauteur du tube. Ces jauges ont des longueurs de 6cm et
des resistances de 120.
Procedure des essais. On applique les charges quasi-statiques en plusieurs cycles. Lamplitude
maximale de la force appliquee est 10kN. La vitesse de deplacement de la barre transversale est
0.5mm/minute et les cycles de chargement sont controles par le logiciel AUTOTRAC. Les signaux
dacquisition sont recuperes en utilisant le logiciel LABVIEW par une voie pour la force appliquee,
une voie pour les capteurs de deplacement, deux voies pour les jauges longitudinales et deux voies
pour les jauges transversales. La figure 6.5 presente une courbe donnant la relation entre la force
appliquee et le deplacement vertical dans les 20 premiers cycles.
Calcul des parametres du ballast Considerons le ballast comme un materiau elastique iso-
trope homogene, les caracteristiques du ballast sont determinees en utilisant la formule de com-
6.4 Identification du comportement statique du ballast 137
1000
-1000
-2000
-3000
-4000
force (N)
-5000
-6000
-7000
-8000
-9000
-10000
-0.0025 -0.002 -0.0015 -0.001 -0.0005 0
deplacement (m)
-1000
-2000
-3000
-4000
force (N)
-5000
-6000
-7000
-8000
-9000
-10000
-0.0024 -0.0023 -0.0022 -0.0021 -0.002 -0.0019 -0.0018 -0.0017 -0.0016
deplacement (m)
portement elastique 3D :
= (tr)1 + 2 (6.75)
L
33 = (6.76)
L
138 6. Modelisation du comportement 3D du ballast
Comme la force appliquee est symetrique et le ballast est suppose homogene, la deformation est
plane dans la section (A-A) et les deformations radiales ( rr ) et tangentes ( ) sont identiques 3 .
Supposons que le deplacement soit continu entre le ballast et le tube, ces deformations sexpriment
a partir de celles du tube PVC quon mesure par deux jauges transversales ( trans ) :
rr = = trans (6.77)
(a) (b)
33 te
rr
rr
rr
33
te
En notant T la force totale de frottement sur la surface interieure du tube, la contrainte verticale
sur le ballast dans la section A-A est :
P T
33 = (6.78)
Aballast
ou Aballast est laire de la section de ballast. La force de frottement T comprime le tube PVC et
provoque une deformation longitudinale quon mesure par deux jauges longitudinales :
PV C
T = Ever longi AP V C (6.79)
on en deduit :
PV C
P Ever PV C
longi A
33 = (6.80)
Aballast
La pression rr va causer une tension t uniforme (en utilisant lhypothese de tube mince) dans
lepaisseur du tube de PVC (t = re rr ) et une deformation trans quon mesure par deux jauges
transversales. Elle est alors determinee :
te E P V C trans e
=rr = (6.81)
r r
On connat toutes les composantes de letat de deformation et de contrainte dans le ballast, en
utilisant (6.75) :
33 = tr + 233 (6.82)
rr = tr + 2rr (6.83)
deux coefficients de Lame sont definis :
33 rr 33 233
= ; = (6.84)
2(33 rr ) tr
Le module dYoung et le coefficient de Poison sexpriment :
(3 + 2) E 2
E= ; = (6.85)
+ 2
On peut donc calculer et a chaque pas du chargement. On suppose que est une fonction de
la dilatation = (I1 ) avec I1 = tr et dans chaque direction est une fonction de la deformation
= (i ).
3
En supposant que la deformation radiale (rr ) dans la section A-A est constante, rr peut sexrimer par :
1 u ur
rr = u
r
= urr . On a : = + = rr comme u = 0 (la deformation est symetrique par rapport a Or).
r r
6.4 Identification du comportement statique du ballast 139
- On constate que le tassement est tres important et diminue de plus en plus apres les cycles de
chargement. On doit donc realiser suffisamment de cycles afin dobtenir un etat stable. En fait,
il faut environ 60-80 cycles de chargement.
- La courbe de force - deplacement au dernier cycle presentee sur la Figure 6.6 a une forme de
type parabolique. La courbure en decharge est plus importante par rapport a la courbe en charge.
Toutefois, la courbe de charge - decharge est fermee et on peut supposer que le ballast est un
materiau elastique non-lineaire.
- Les figures 6.8 et 6.9 tracent les courbes de et en fonction de I 1 et . On constate que la
rigidite du ballast saccroit en fonction de la deformation. En plus, ces courbes I 1 et
sont assez regulieres et ont aussi des formes paraboliques.
100
essai 1
essai 2
90 essai 3
essai 4
valeur moyenne
80
70
lambda (MPa)
60
50
40
30
20
10
-0.004 -0.0035 -0.003 -0.0025 -0.002 -0.0015 -0.001 -0.0005
I1
50
essai 1
essai 2
45 essai 3
essai 4
valeur moyenne
40
35
mu (MPa)
30
25
20
15
10
-0.004 -0.0035 -0.003 -0.0025 -0.002 -0.0015 -0.001 -0.0005
epsilon
- Les valeurs calculees sont tres variees suivant les essais meme avec des conditions aux limites
identiques (le tube et le chargement), on peut remarquer que le caractere mecanique du ballast
depend beaucoup de letat initial des cailloux dans le tube.
- Afin de determiner les formules de et , on prend les valeurs moyennes des essais realises.
On remarque que ces valeurs peuvent etre ajustees avec des polynomes du 4 e ordre dont les
coefficients sont les suivants :
= 7.1 1011 I14 + 3.18 109 I13 + 7.5 106 I12 + 1.54 104 I1 + 18.5 (6.86)
11 4 8 3 6 2 2
2 = 1.8 10 + 9.51 10 + 2.89 10 5.94 10 + 14.2 (6.87)
p = 88.42 103 M P a
d d = 20cm
h = 10cm
p
h
45o
On utilise le modele elastique non-lineaire pour le ballast. La structure est calculee dans les
deux cas delasticite normale et non-tension. Les coefficients et sont pris a partir des resultats
experimentaux (eqs. 6.86,6.87).
Les figures suivantes presentent les resultats calcules dans les cas des lois de comportement
elastique non-lineaire normale (Figures 6.11 - 6.14) et unilaterale (Figures 6.15 - 6.18). Comme la
structure est symetrique, il suffit de visualiser les resultats dans une section verticale. La compa-
raison en chiffre est presentee dans le tableau 6.5.
On constate dans le contour des contraintes majeures (6.11) du cas elastique quil y a une zone
de contraintes positives a cote de la surface laterale. La valeur maximale des contraintes positives
( 1.5 103 M P a) est assez importante (environ 14%) par rapport a la valeur minimale des
contraintes negatives ( 1.0 102 M P a Fig. 6.12). Au contraire, avec le modele non-tension,
il ny a pas du tout de contraintes positives (Fig. 6.15 et 6.16). Les contraintes dans ce modele
sarrangent et la valeur minimale des contraintes negatives est environ 14% plus importante par
rapport au cas elastique.
La structure dans le cas non-tension est evidemment plus souple que lautre. On a une valeur
maximale de 0.0634mm pour le deplacement vertical dans le cas non-tension (Fig. 6.18). Par
6.5 Validation numerique 141
comparaison avec le cas elastique ou elle vaut 0.0414mm, elle est 53% plus grande. Aussi dans
les deplacements lateraux, la difference est de 210%.
Tableau 6.2.
Valeurs maximales de Modele classique Modele non-tension Difference
Contrainte principale positive 1.5 10 3 MPa 0
Contrainte principale negative 1.03 10 2 MPa 1.18 102 MPa 15%
Deplacement vertical 4.14 10 2 mm 6.344 102 mm 53%
Deplacement horizontal 1.38 10 2 mm 4.32 103 mm 210%
Contrainte majeure S1
.00485642 .00133286
.00415171 .00062814
.003447 7.6568E5
.00274228 .00078128
.00203757 .00148599
.00133286 .00219071
Contrainte mineure S3
.0103319 .00484657
.00923486 .00374949
.00813779 .00265242
.00704072 .00155535
.00594364 .00045827
.00484657 .0006388
Deplacement U
0. 6.8872E6
1.3774E6 8.2647E6
2.7549E6 9.6421E6
4.1323E6 1.102E5
5.5098E6 1.2397E5
6.8872E6 1.3774E5
Deplacement W
4.141E5 2.0705E5
3.7269E5 1.6564E5
3.3128E5 1.2423E5
2.8987E5 8.2821E6
2.4846E5 4.141E6
2.0705E5 0.
.00368012 .00173384
.00329086 .00134458
.00290161 .00095533
.00251235 .00056607
.0021231 .00017682
.00173384 .00021244
.0118205 .00583958
.0106244 .00464339
.00942816 .00344719
.00823197 .002251
.00703578 .00105481
.00583958 .00014139
0. 2.3994E5
4.7988E6 2.8793E5
9.5976E6 3.3592E5
1.4396E5 3.8391E5
1.9195E5 4.3189E5
2.3994E5 4.7988E5
Deplacement W
6.3468E5 2.9246E5
5.6624E5 2.2401E5
4.9779E5 1.5557E5
4.2935E5 8.7126E6
3.609E5 1.8681E6
2.9246E5 4.9763E6
6.6 Conclusion
On a propose un modele continu macroscopique des materiaux ballast. Vu les resultats dessai
obtenus, on peut dire que le comportement du ballast est du type elastique non-lineaire en com-
pression. En revanche, la loi de comportement est introduite pour que les deformations en tension
dans certaines directions ne causent aucune contrainte normale dans cette direction. Une fonction
denergie de deformation est proposee en modifiant celle du cas elastique classique. La resolution
numerique est faite par une procedure classique delements finis pour le probleme elastique non-
lineaire qui est ecrit dans la base principale des deformations.
Cette idee peut etre utilisee pour les autres materiaux de type granulaire (comme le sable
...). En plus, avec la fonction denergie de deformation proposee, on peut facilement considerer les
cas plus compliques des materiaux avec de la visco-elasticite unilaterale ou de lelasto-plasticite
unilaterale.
144 6. Modelisation du comportement 3D du ballast
Chapitre 7
7.1 Introduction
On a vu dans le chapitre 2 que le probleme dun massif multi-couche soumis a une charge mobile
peut etre resolu par une methode semi-analytique. Cette methode donne des resultats exacts mais
nest valable que dans le cas elastique. Dans les cas ou on a un comportement non-lineaire des
materiaux (qui nest pas du tout negligeable pour des materiaux comme le ballast), la methode
des elements finis doit etre envisagee. Le probleme se ramene a resoudre un systeme dequations
differentielles en temps avec les conditions initiales :
ou M et K sont respectivement les matrices de masse et de rigidite, F est le vecteur des forces
exterieures qui nest pas fixe dans lespace. Un tel vecteur force demande un gros maillage pour
pouvoir saffranchir du regime transitoire ce qui est presque impossible dans le cas tridimensionnel.
Les parties suivantes presentent une methode des elements finis modifiee qui permet de passer
directement dans le regime permanent et de traiter statiquement le probleme dans un domaine
mobile avec une modification de la matrice de rigidite.
Le probleme consiste a determiner quels sont les champs mecaniques (deplacement, vitesse,
acceleration) stationnaires dans la structure. On note u(x, t) le vecteur deplacement : (x, t)
146 7. Massif 3D non-lineaire soumis a une charge mobile.
]0, T [ ou T est lintervalle de temps du calcul. On peut ecrire lequation dequilibre dans et
les conditions aux limites :
2 u(x, t)
Divx (x, t) = 0 x
t2
(7.4)
(x, t) n = f (x, t) x
u(x, t) = 0 x u
ou est la masse volumique, n est le vecteur normal a la surface , designe le tenseur des
contraintes qui est relie au tenseur des deformations par une loi de comportement qui, dans le
cas elastique, secrit ij () = Cijkl kl avec (x, t) = 12 [u(x, t) + t u(x, t)].
. f
/
repere mobile
z
0 (x vt)(y)(z)
repere fixe y
z x
O
y
O x vt
u : 0 0
=
Figure 7.1.
En remplacant ces relations dans (7.4), le systeme dequations dans le domaine mobile secrit :
2 u (x , t) 2 u (x , t) 2
2 u (x , t)
2v + v Divx (x , t) = 0 x
t2 tx x 2
(7.8)
(x ) n = f 0 (x )(y)(z) x
u (x ) = 0 x u
Le probleme consiste donc a resoudre un systeme dequations aux derivees partielles du 2 e ordre
par la methode des elements finis. Afin de simplifier la notation, on va supprimer les dans les
expressions qui suivent.
7.2 Equations du probleme 147
Nous utilisons la procedure classique des elements finis pour resoudre le systeme (7.8). Lex-
pression de la formulation faible de ce probleme secrit :
Z
2 u(x, t) 2 u(x, t) 2
2 u(x, t)
W = < 2v + v Divx (x), u(x) > d = 0
t2 tx x2
u(x) C ad (7.9)
ou <, > represente le produit scalaire et u est la fonction test qui est definie dans lespace des
deplacements admissibles C ad .
2
Pour simplifier la notation, on note u u u
t = u, t2 = u et x = u,x . En utilisant la formule de
Green dans lintegration (7.9) et en tenant compte des conditions aux limites, on obtient :
Z Z Z
W = u u d 2v u,x u d v 2 u,x u,x d
Z Z (7.10)
+ ( : u) d f u d = 0
X
Le domaine est represente par un ensemble de sous domaines = e et lintegration (7.10)
P P R R e
est definie par : W = e W = e [ e (.) + e (.)]. Lintegration dans chaque element e peut etre
evaluee en fonction des deplacements aux noeuds U e par lapproximation :
n
X
e
u (x, t) N i (x)U ei (t) = N (x)U e (t) (7.11)
i=1
ou N est la fonction dinterpolation. La fonction test du type de Galerkin utilise la meme fonction
dinterpolation que ue : ue = N U e .
Rappelons que pour la mise en oeuvre de la methode des elements finis, le produit ( : ) est
habituellement represente par t qui peut sexprimer en fonction de U grace a (7.11) [51] :
= {11 22 33 12 13 23 }t = D (7.12)
t
= {11 22 33 212 213 223 } = BU (7.13)
Lassemblage des sous domaines e (7.14) mene au systeme dequations dynamique suivant :
ou :
XZ
M est la matrice de masse : M = N t N d.
e e
148 7. Massif 3D non-lineaire soumis a une charge mobile.
XZ
C est la matrice damortissement mobile : C = 2v N t N ,x d.
e e
Remarque 7.2.1
1. On peut montrer facilement que [Kv ] et [Ks ] sont deux matrices symetriques et positives. La symetrie
de [K] est donc toujours assuree.
2. La matrice damortissement C est non-symmetrique.
3. Lequation (7.15) presente tres bien linfluence de la vitesse de la charge mobile sur la structure. Le
probleme mobile a une rigidite plus faible que dans le cas statique. Dans le cas statique, le vecteur
deplacement est obtenu simplement par {U} = [K]1 {F}. Plus la vitesse est elevee, plus la rigidite
mobile est faible, et plus le champ de deplacement est important.
4. La positivite de [K] depend de la valeur de la vitesse v. La matrice de rigidite [K] pourrait ne pas etre
positive avec une vitesse v suffisamment grande. Dans ce cas, la forme faible (7.9) nest plus valable
(theoreme de Lax-Migram).
Dans le cas ou lamplitude de la charge est constante, la solution stationnaire dans le repere
mobile devient statique, et le systeme dequations delements finis se reduit a :
K U = F (7.16)
Dans le cas nonlineaire, une procedure iterative doit etre envisagee. On peut utiliser la methode
de Newton-Raphson classique pour ce probleme de charge mobile (7.16). Supposons quon a une loi
de comportement qui secrit sous la forme : {} = [D()]{}, la procedure de calcul est resumee
dans le schema suivant :
(a) A literation i : U (i) et (i) :
Calcul du vecteur residu :
Z Z
R(i) = B t (i) d v 2 N t,x N ,x U (i) d F (7.17)
ou tol est une tolerance autorisee. Si la condition (7.21) est verifiee, la solution finale est
{U } = {U }(i+1) . Sinon, on revient en (a) et on calcule la nouvelle iteration.
On voit bien que par rapport au cas statique, la vitesse nintervient dans la matrice tangente de
rigidite que comme une constante. Le calcul du comportement des materiaux peut donc etre realise
de facon independante.
La matrice de rigidite modifiee dans le domaine mobile a ete implementee dans le code delements
finis CESAR du LCPC.
7.3 Validation
On presente une validation numerique dans le cas ou il y a une force qui se deplace sur la
surface libre dun demi espace elastique infini. Comme ce probleme est symetrique par rapport
au plan Ox, Oz, on peut enlever une moitie du modele (la partie y < 0) et ajouter la condition
uy = 0. La structure est modelisee dans un domaine de dimension (16, 16) (0, 16) (16, 0) par
un maillage de 9581 noeuds et de 2000 elements cubiques de 20 noeuds. Les noeuds aux surfaces
exterieures sont fixes. Les parametres du materiau utilise sont : module dYoung E = 130M P a,
coefficient de Poisson = 0.3. Une force ponctuelle de 1000N est appliquee au point (0, 0, 0).
z
y
0
x
uy = 0
Ci-dessous sont les solutions des deplacements verticaux au niveau de la profondeur de 0.35m
dans le plan y = 0. La figure (7.3) presente une comparaison entre la solution analytique et la
solution par elements finis dans le cas de la vitesse v = 50m/s. On obtient une petite difference
entre ces deux solutions qui vient du fait quon a utilise une condition fixe aux limites.
150 7. Massif 3D non-lineaire soumis a une charge mobile.
1.5e-06
-2e-06
1e-06
-4e-06
5e-07
-6e-06
u1 (m)
u3 (m)
0
-8e-06
-5e-07
-1e-05
-1e-06
-1.2e-05
-1.5e-06
MEF MEF
analytique analytique
-2e-06 -1.4e-05
-4 -2 0 2 4 -4 -2 0 2 4
X (m) X (m)
Dans la figure (7.4), on compare les deplacements verticaux dans les cas ou la force se deplace
avec les vitesses de 50m/s, 100m/s,150m/s,180m/s. Cette comparaison montre bien leffet tres
important de la vitesse de la charge sur la structure.
On remarque aussi que le probleme ne peut pas etre resolu pour les vitesses plus elevees comme
pour ce massif, la vitesse de Rayleigh vaut 190 m.s 1 .
5e-06 5e-06
4e-06
0
3e-06
-5e-06
2e-06
Deplacement ux (m)
Deplacement uz (m)
-1e-05
1e-06
0 -1.5e-05
-1e-06
-2e-05
-2e-06
-2.5e-05
-3e-06
v = 50m/s v = 50m/s
v = 100m/s -3e-05 v = 100m/s
-4e-06
v = 150m/s v = 150m/s
v = 180m/s v = 180m/s
-5e-06 -3.5e-05
-4 -2 0 2 4 -4 -2 0 2 4
X (m) X (m)
On sinteresse maintenant au probleme dun massif non-lineaire soumis a une charge mobile.
Pour ce faire, on va introduire une non linearite de type unilateral (quon a presente dans le chapitre
precedent) dans la structure consideree.
On suppose que le materiau est elastique lineaire unilateral avec un module dYoung en tension
(Et ) qui vaut la moitie de celui en compression (E c = 130 M P a).
La figure 7.5 presente les resultats de deplacement horizontal et vertical en comparant avec ceux
des cas lineaires : soit avec une rigidite E t , soit avec une rigidite Ec . On a obtenu une solution de
deplacement vertical qui se trouve entre les deux solutions lineaires, puisque ce massif non-lineaire
est plus souple que celui avec Ec et est plus rigide que lautre.
7.4 Couplage dune poutre et dun demi-espace a deux couches 151
2e-06 -1e-05
1e-06
-1.5e-05
u3 (m)
u1 (m)
0
-2e-05
-1e-06
-2e-06 -2.5e-05
-3e-06
unilateral -3e-05 unilateral
-4e-06 E = 130 MPa E = 130 MPa
E = 65 MPa E = 65 MPa
-5e-06 -3.5e-05
-4 -2 0 2 4 -3 -2 -1 0 1 2 3
X (m) X (m)
On realise dans cette partie des calculs sur un probleme plus realiste. Le modele de calcul est
un couplage entre une poutre, qui represente le rail, et un demi-espace constitue de deux couches,
qui representent la couche de ballast et sol.
7.4.1 Maillage
Le maillage de calcul est presente dans la figure 7.6. Comme le probleme est symetrique par
rapport au plan Oxz, on ne maille que la partie y 0. Le maillage se compose de trois groupes
delements de differents materiaux pour le rail, la couche de ballast et le sol.
M
ZOO
Le rail. Pour simplifier la geometrie de la section-I du rail mais garder la meme repartition de
charge sur le ballast, on introduit une section de forme rectangulaire qui a toujours la meme
largeur et aussi le meme moment dinertie I. Appelons la largeur et la hauteurqde cette section
3 12I
rectangulaire par b et h, elles sont determinees en fonction de S et I par : h = b et la masse
volumique utilisee pour le rail est egale a 0 = Sbh
rail
. Dans ce calcul ou le rail quon utilise a
4
b = 15 cm, S = 60.34 kg/m et I = 3060 cm , on peut utiliser une section rectangulaire dont
h = 13.5 cm, b = 15 cm et 0 = 2980 kg/m3 (1 ).
Dans CESAR, cette poutre est maillee par des elements volumiques de type cubique de 20 noeuds
et on utilise quatre elements pour chaque section. Dans le calcul, on suppose quelle reste toujours
elastique lineaire.
Le ballast. La couche de ballast est modelisee comme une couche infinie dont lepaisseur est de
30 cm. Cest dans cette couche que la non linearite sera introduite. Cette couche est modelisee
par des elements cubiques de 20 noeuds.
Pour les calculs lineaires, le module dYoung du ballast est egal a 130 M P a. Pour les calculs
non-lineaires, on utilise les polynomes de et quon a identifies dans le chapitre precedent.
elastique lineaire elastique lineaire unilateral elastique nonlineaire elastique lineaire unilateral
Le sol. Le sol est aussi suppose etre elastique lineaire et est maille par des elements prismes de
15 noeuds. Le module dYoung du sol est pris egal a E s = 100 M P a.
Le maillage a 16761 noeuds au total avec 4088 elements (2800 elements cubiques de 20 noeuds et
1288 elements prismes de 15 noeuds).
La charge est une force constante damplitude 170 kN appliquee sur la poutre.
Les resultats sont obtenus en realisant des calculs avec les quatre cas de comportement du
ballast : elastique lineaire et elastique nonlineaire, avec ou sans comportement unilateral.
On visualise dans la figure 7.7 les contours des contraintes principales maximales dans le plan
y = 0. On peut constater que la zone ou il existe des contraintes principales positives dans le
modele lineaire va disparatre en utilisant la loi unilaterale. Par contre, les contraintes negatives
sont plus importantes dans le modele unilateral.
La figure 7.8 presente les comparaisons des deplacements verticaux au niveau z = 15 cm (la
moitie de lepaisseur de la couche de ballast) dans la surface y = 0. On constate que les deplacements
verticaux avec le comportement unilateral sont plus importants, surtout dans le cas non-lineaire.
Dans le cas non-lineaire, le deplacement statique du a la gravite est beaucoup plus important avec
le modele unilateral.
1
On peut aussi utiliser une section rectangulaire equivalente en gardant le meme , mais elle aurra une largeur
tres petite qui ne peut pas bien representer la repartion de la charge
7.4 Couplage dune poutre et dun demi-espace a deux couches 153
Contrainte majeure S1
148091. 35870.8
125647. 13426.8
103203. 9017.25
80758.9 31461.3
58314.9 53905.3
35870.8 76349.4
Contrainte majeure S1
166695. 54519.1
144260. 32084.
121825. 9648.81
99389.4 12786.3
76954.3 35221.5
54519.1 57656.6
v = 50 m/s v = 50 m/s
-0.7 -2
-0.8
-3
-0.9
-4
-1
-1.1 -5
u3 (mm)
u3 (mm)
-1.2 -6
-1.3
-7
-1.4
-8
-1.5 lineaire non-lineaire
lineaire unilateral non-lineaire unilateral
-1.6 -9
-10 -5 0 5 10 -10 -5 0 5 10
X (m) X (m)
-3 -4
-3.2 -4.5
-3.4 -5
-3.6 -5.5
u3 (mm)
u3 (mm)
-3.8 -6
-4 -6.5
-4.2 -7
-4.4 -7.5
v = 0 m/s v = 0 m/s
-4.6 v = 25 m/s -8 v = 25 m/s
v = 50 m/s v = 50 m/s
-4.8 -8.5
-4 -2 0 2 4 -4 -2 0 2 4
X (m) X (m)
Avec les deplacements obtenus, on peut aussi determiner les accelerations par la relation u(x) =
v 2 ,xx u(x). La figure 7.10 presente des comparaisons des accelerations calculees avec les modeles
differents. Les derivees en x sont calculees directement a partir des vecteurs de deplacement en
utilisant des outils numeriques. Les accelerations obtenues avec les modeles lineaires sont nettement
inferieures a celles calculees avec les modeles non-lineaires. Cette remarque se signifie, quavec un
v = 50 m/s v = 50 m/s
6 20
lineaire nonlineaire
5 lineaire unilateral nonlineaire unilateral
15
4
acceleration verticale (ms2)
acceleration verticale (ms )
2
3 10
2
5
1
0 0
1
5
2
3 10
10 5 0 5 10 10 5 0 5 10
X (m) X (m)
ballast bien bourre dont le comportement est assez rigide et lineaire, lacceleration verticale est
moins importante.
7.5 Conclusion
On a presente une formulation par elements finis pour resoudre les problemes de charge mobile
en regime permanent par un changement de coordonnees. Lavantage de cette methode est que lon
peut traiter le probleme avec une charge fixe comme en statique. Cette methode nest valable que
dans les cas ou la vitesse de la charge reste subsonique.
Le modele nest pas encore etudie en detail, mais on peut faire les remarques suivantes :
Dans le cas elastique, cette methode est valable seulement dans les cas ou la vitesse de la charge
est inferieure a une vitesse critique qui correspond a la plus petite vitesse donde de Rayleigh
dans le milieu.
Dans le cas unilateral, la vitesse critique diminue. On remarque aussi que meme si la rigidite en
tension dans le ballast est nulle, il y a toujours des ondes de cisaillement qui peuvent se propager
de facon normale.
156 7. Massif 3D non-lineaire soumis a une charge mobile.
Conclusions et Perspectives
Conclusions
Les problemes non-lineaires sont traites en utilisant la methode des elements finis, soit dans
le repere fixe (a), soit dans le repere mobile (b). On a vu que le cout des calculs stationnaires
dans le repere mobile est plus faible que les calculs transitoires dans le repere fixe. Deux approches
differentes pour le calcul dans le repere mobile sont utilisees. La premiere (b1) transforme directe-
ment les equations dans le repere mobile par un changement de variable. Elle est plus simple pour
la mise en oeuvre numerique mais nest valable que dans les cas subsoniques. La deuxieme (b2)
realise un changement de variable a partir des equations discretisees. Elle permet de resoudre le
probleme avec toutes les vitesses mais a ete validee uniquement pour les modeles 1D. On presente
dans le tableau ci-dessous les problemes et les approches correspondantes possibles, ainsi que celles
utilisees dans ce memoire.
Par ailleurs, dans le calcul des modeles 1D de poutres, on a introduit une condition aux limites
absorbante qui peut amortir les ondes en utilisant une couche dont lamortissement crot en fonction
de la distance au point considere et a la source.
Perspectives
Le modele semi-analytique est tres utile pour des etudes preliminaires grace a sa rapidite en
temps de calcul meme pour une structure 3D. Pour lapplication aux voies ferrees, une etude plus
precise peut etre developpee en ajoutant dans le probleme presente dans le chapitre 2 une poutre
dEuler-Bernouilli (ou de Timoshenko). Les appuis periodiques (les blochets) peuvent eventuellement
etre introduits afin davoir un modele lineaire complet. Autrement, cette methode peut etre aussi
utilisee pour des modeles viscoelastiques de routes. Un travail possible par la suite est dintroduire
la loi de Huet-Sayleigh pour decrire le comportement visco-elastique des enrobes et de lintegrer
dans le logiciel ALIZE.
Le probleme non-lineaire doit etre etudie plus en detail pour les massifs. A partir des resultats
obtenus dans cette these, on peut developper un modele numerique sur CESAR incluant les com-
portements non-lineaires pour le ballast, les semelles et les sols, dans un premier temps dans le
cas bidimensionnel. Le premier modele est donc un rail pose sur des semelles visco-elastiques non-
lineaires elles-memes reposant sur des traverses continues. Le tout est pose sur du ballast et sur un
sol multicouche. Ensuite, les traverses continues peuvent etre remplacees par des appuis discrets
pour connaitre linfluence des effets des appuis periodiques dans le cas dune charge mobile. Ces
modeles permettent aussi de savoir sil est possible de remplacer le comportement multicouche bi-
dimensionnel du sol par un comportement unidimensionnel equivalent. Dans le cas dune reponse
positive, on fournira un moyen de determiner les parametres 1D equivalents aux comportement
2D. Les memes problematiques peuvent etre etudiees pour le modele tridimensionnel.
En parallele, on peut effectuer des mesures sur le banc dessais dabord pour obtenir la deflexion
sous chargement statique qui sera comparee aux resultats des modeles 1D et 2D. En dynamique,
on peut faire des mesures de fonctions de transfert entre un signal dexcitation et le signal mesure
par des accelerometres en differents points de la structure. Enfin on peut solliciter la structure
avec le signal M applique habituellement pour les essais de fatigue et on peut mesurer de nouveau
lacceleration en differents points de la structure. Des calculs de structures avec les differents
modeles pour ces trois types dexcitations peuvent etre compares avec les resultats de mesures.
Conclusions & Perspectives 159
On a etudie ici linfluence des non-linearites dans les lois de comportement en restant cependant
dans le domaine des petites deformations. Il conviendrait dans une etude ulterieure de traiter le
cas des grandes deformations.
160 Conclusions & Perspectives
Annexes
Annexe A
avec :
x2n x0
h = (A.2)
2n
Xn
C2n = f2i cos(tx2i ) 21 [f2n cos(tx2n ) + f0 cos(tx0 )]
i=0
n
X
C2n1 = f2i1 cos(tx2i1 ) (A.3)
i=1
Xn
0 (3)
S2n1 = f2i1 sin(tx2i1 ) (A.4)
i=1
1 sin(2) 2 sin2
() = + (A.5)
2 2 3
1 + cos2 sin(2)
() = 2 (A.6)
2 3
sin cos
() = 4 2 , (A.7)
3
On presente ici quelques resultats analytiques du probleme dun demi-espace soumis aux charges
statiques a sa surface libre (ref. [56, chap. 3]).
r x
A0
y z
Figure B.1.
p
= x2 + y 2 + z 2 (B.1)
p
r = x2 + y 2 (B.2)
Les deplacements et les contraintes a un certain point A(x, y, z) dans le demi-espace sont donnes
par :
166 B. Charge statique sur un demi-espace elastique
P xz x
ux = (1 2) (B.3)
4G 3 ( + z)
P yz x
uy = (1 2) (B.4)
4G 3 ( + z)
2
P z 2(1 )
uz = + (B.5)
4G 3
P 1 2 z x2 y 2 zy 2 3zx2
xx = 1 + 3 5 (B.6)
2 r2 r2
2 2 2
P 1 2 z y x zx 3zy 2
yy = 1 + 3 5 (B.7)
2 r2 r2
3P z 3
zz = (B.8)
2 5
P 1 2 z xy xyz 3xyz
xy = 1 3 5 (B.9)
2 r2 r2
3P xz 2
xz = (B.10)
2 5
3P yz 2
yz = (B.11)
2 5
4(1 2 )pa
uz = E(r/a) si r<a (B.12)
E
4(1 2 )pa
uz = [E(a/r) (1 a2 /r 2 )K(a/r)] si r>a (B.13)
E
ou on note K et E les integrales elliptiques completes de premier et second type, respectivement.
Au centre du cercle, r = 0, E(0) = 2 et a cote du cercle r = a, E(1) = 1, donc :
2(1 2 )pa
(uz )0 = (B.14)
E
4(1 2 )pa
(uz )a = (B.15)
E
Le deplacement moyen du cercle de charge est :
16(1 2 )pa
u = (B.16)
3E
Annexe C
Matrices elementaires
13 1 1
11 h 11 1 1
11 h
2 7 1 1
12 h
0
9 1 1
11 h
0
3 1 1
12 h
0 0 0 0 0
35
1
11 111 h2 1 210
1
11 h
3 1 20
1
12 h
2
0
13 1 70
1
11 h
2
1 1
11 h
3 1 20
1
12 h
2
0 0 0 0 0
1
210
1 105
1 20
1 1 420
1 140
1 30
1
7 121 h 1
21 h
2 1
22 h
1
23 h 3 1
21 h
1
21 h
2 1
22 h
1
23 h
0 0 0 0
20 20
1 3
1 3 20 30
1 6
1 6
1 1 1 1
32 h 33 h 32 h 33 h
0 0 0 0 0 0 0 0
Mea =
1 1 1 3 3
1 1 1 6 6
9 111 h 13 1
11 h
2 3 1
12 h
0
13 1
11 h
11 1
11 h
2 7 1
12 h
0 0 0 0 0
13701 h2
11 1
1 420
1
11 h
3
1 20
1
12 h
2
0
35
11 111 h2
210
1 1 1
11 h
210
3
1 20
1
12 h
2
0 0 0 0 0
1
420
3 1h
21 1 1
140
21 h
2 1 1
30
22 h 1 1
23 h 1
7 121 h
1 105
1
21 h
2 1 1
20
22 h 1 1
23 h
0 0 0 0
1 1 1 1
20 30 6 6 20 20 3 3
1 1 1 1
32 h 33 h 32 h 33 h
0 0 6 6 0 0 3 3 0 0 0 0
(C.1)
12 6
h3 h2
0 0 h123 6
h2
0 0 0 0 0 0
6 4
h2 h 0 0 h62 2
h 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
Keaf = A211
12 (C.2)
h3 h62 0 0 12
h3
h62 0 0 0 0 0 0
6 2
0 0 h62 4
0 0 0 0 0 0
h2 h h
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0
13 2 1
11 h 11 2 1
11 h
2 71
12 h 2 0
9 2 1
11 h
0
31
12 h 2 0
0 0 0 0
11 2 35
1
11 h
2 2 210
1
11 h
3 1
2
20
12 h
2
0
13 2 70
1
11 h
2
2 1
11 h
3 1
2
20
12 h
2
0
2 210
2 105 20
2 2 2 420
2 140 30
2 2
7 1 1 2 1 1
3 1 1 2 1 1
21 h 21 h 22 h 23 h 21 h 21 h 22 h 23 h
0 0 0 0
20 20 3
2 2 3 20 30 6
2 2 6
1 1 1 1
32 h 33 h 32 h 33 h
0 0 0 0 0 0 0 0
Meb =
2 2 3
2 3
2 2 6
2 6
9 1
13 1 2 3 112 h 13 1
11 1 2 7 112 h
11 h 11 h 11 h 11 h
0 0 0 0 0 0
70
2 2 420 20
2 35
2 2 210 20
2
13 111 h2 1 3 1
h2 11 111 h2 1 3 1
h2
11 h 11 h
0 0 0 0 420 12 0 210 12 0
2 2 2 2 2 2 2 2
140 30 105 20
3 121 h 1 2 1 1
7 121 h 1 2 1 1
21 h 22 h 23 h 21 h 22 h 23 h
0 0 0 0
2 2 2 2
20 30 6 6 20 20 3 3
1 1 1 1
32 h 33 h 32 h 33 h
0 0 0 0 0 0 6 6 0 0 3 3
(C.3)
168 C. Matrices elementaires
0 0 0 0
2 2
11 2 2
11 h
2 2
12
0
2 2
11
2 2
11 h 2 2
12
0
0 0 0 0
2 2
2
11 h
10
0
2 2
2
12 h
0
2 2
2
11 h
2 10
2
11 h
2 2 2
2
12 h
0
2 10
2 2 12
2 2 10
2 60
2 12
2
2 2 2 2 2 2 2 2
21 h 21 h
0 0 0 0 21
22
23 21
22 23
2 12
2 2
2 2 2 12
2 2
2 2
2 2 2 2
0 0 0 0 0 0 32
33
0 0 32 33
e
2 2 2 2 2
2 2 2 2 2
Cb = 2 2 2 2 2 2 (C.4)
11 h 11 h
0 0 0 0 11
12
0 11 12
0
2 2
2 10
2 2
2 2 10
2 2
2 2 2 2 2 2
11 h 11 h 12 h 11 h 12 h
0 0 0 0 0 0 0
2 2 2 2 2 2 2 2
10 60 12 10 12
2 2 2 2 2 2 2 2
21 h 21 h
0 0 0 0 21
22
23 21 22 23
2 2 2 2
2 12 2 2 2 12 2 2
2 2 2 2
0 0 0 0 0 0 2
32
2
33
0 0 2
32
2
33
0 0 0 0
6 2 3
11 2 3
11 2 3
12
0
6 2 3
11 2 3
11
2 3
12
0
0 0 0 0
2 5h
3
11 2 2 10
3
11 h
h
0 0
2 5h
3
11
2 10
3
11 h
0
h
0
2 10 15
2 2 2 10 30
2 2
3 3 3 3 3 3
0 0 0 0 21
0 22 23
21
0 22
23
h
2 h
2 h h
2 h
2 h
3 3 3 3
0 0 0 0 0 0 32 32
0 0 32
32
Kebv
2 2 2 h h
2 2 2 h h
= 3 3 3 3 3 3 (C.5)
6 6
0 0 0 0 11
11
12
0 11
11 12
0
2 5h
2 10 h
2 5h
2 10 h
3 3 3
2 3
11 h 11 h
0 0 0 0 11
0 0 11
0 0
2 2 2 2 2 2
10 30 10 15
3 3 3 3 3 3
0 0 0 0 21
0 22
23 21
0 22 23
2 2 2 2
h h h h h h
3 3 3 3
0 0 0 0 0 0 h
32
h
32
0 0 h
32
h
32
12 6
0 0 0 0 h3 h2
0 0 h123 6
h2
0 0
6 4
0 0 0 0 h2 h 0 0 h62 2
h 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
Kebf = B 411
(C.6)
0 0 0 0 h123 h62 0 0 12
h3
h62 0 0
6 2
0 0 0 0 h2 h 0 0 h62 4
h 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
Bibliographie
[1] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions. Dover Publication,
Inc., 9th edition, 1972.
[2] J.D. Achenbach. Wave Propagation in Elastic Solids. North-Holland, 1973.
[3] L. Achimastos. Appreciation de letat structurel dune chaussee souple a partir des observa-
tions de degradations de surface. PhD thesis, Ecole Nationale des Ponts et Chaussees, March
1998.
[4] A. Alaoui and T. Naciri. Les voies ballastees. Rapport EUROBALT. CERAM, 1995.
[5] G. Alfano, L. Rosati, and N. Valoroso. A numerical strategy for finite element analysis of no-
tension materials. International Journal for Numerical Methods in Engineering., 48 :317350,
2000.
[6] J. Alias. La voie ferree - Techniques de construction et dentretien. Eyrolles, deuxieme
edition, 1984.
[7] L. Andersen, S.R.K. Nielsen, and P.H. Kirkegaard. Finite element modelling of infinite Euler
beams on Kelvin foundations exposed to moving loads in convected co-ordinates. Journal of
Sound and Vibration, 241(4) :587604, 2001.
[8] D. Aubry, D. Clouteau, and A. Modaressi. Interaction dynamique sol-structure. In Calcul
tridimensionnel en geotechnique. Presse des Ponts et Chaussees, 1998.
[9] E. Azanka. Ecoulements granulaires bidimensionnels sur un plan incline. PhD thesis, Ecole
Nationale des Ponts et Chaussees., January 1998.
[10] K. Bagi. Stress and strain in granular assemblies. Mechanics of Materials, 22(165177), 1996.
[11] M.C.M. Bakker, M.D. Verweij, B.J. Kooij, and H.A. Dieterman. The traveling point revisited.
Wave Motion, 29 :119135, 1999.
[12] M. Balsan. Pour un modele mathematique de la voie ferree moderne. PhD thesis, Ecole
Nationale des Ponts et Chaussees., 1980.
[13] M. Bats-Villard. Influence des defaults de liaison sur le dimensionnement et le comportement
de chaussees. PhD thesis, Universite de Nantes, 1991.
[14] P.M. Belotserkovskiy. On the oscillations of infinite periodic beams subjected to moving
concentrated force. Journal of Sound and Vibration, 193(3) :705712, 1996.
[15] P.M. Belotserkovskiy. Forced oscillations of infinite periodic structures. Applications to rail-
way track dynamics. Vehicle System Dynamics Supplement, 28 :85103, 1998.
[16] V. Bodin. Comportement du ballast des voies ferrees soumises a un chargement vertical et
lateral. PhD thesis, Ecole Nationale des Ponts et Chaussees, June 2001.
[17] V. Bodin, P. Tamagny, K. Sab, and P.-E. Gautier. Frequence critique de sollicitation lors du
tassement a grand nombre de cycles dun milieu granulaire. Colloque Physique et Mecanique
des Materiaux Granulaires. Champs-sur-Marne, 5-7 Septembre, 2000. Actes des journees
scientifiques du LCPC, tome 2, pages 291296, 2000.
170 BIBLIOGRAPHIE
[38] E.S. Eringen and A.C. Suhubi. Elastodynamics. Vol. 2 - Linear theory. Academic Press,
1975.
[39] L. Fryba. Vibration of Solids and Structures under Moving Loads. Thomas Telford, 3rd
edition, 1999.
[40] S. Garvilov. Non-stationary problems in dynamics of a string on an elastic foundation sub-
jected to a moving load. Journal of Sound and Vibration, 222(3) :345361, 1999.
[41] M. Geradin and D. Rixen. Theorie des vibrations. Application a la dynamique des structures.
Masson, deuxieme edition, 1996.
[42] J. Ghabussi and R. Barbosa. Three-dimensional discrete element method for granular ma-
terials. J. Num. & Anal. Methods Geomech., 14 :451472, 1990.
[43] V.T. Granik and M. Ferrari. Microstructural mechanics of granular materials. Mechanics of
Materials, 15(301322), 1993.
[44] H. Grundmann, M. Lieb, and E. Trommer. The response of a layered half-space to traffic
loads moving on its surface. Arch. Appl. Mech., 69 :5567, 1999.
[45] L. Gry. Modelisation du comportement dynamique dune voie TGV pour la reduction du bruit
de roulement. PhD thesis, Ecole Centrale Paris., November 1995.
[46] N. Guerin. Approche experimentale et numerique du comportement du ballast des voies
ferrees. PhD thesis, Ecole Nationale des Ponts et Chaussees., November 1996.
[47] J.-F. Hamet. Railway noise : Use of Timoshenko model in rail vibration studies. Acta
Acustica, 85 :5462, 1999.
[48] M.S.A. Hardy. The generation of waves in infinite structures by moving harmonic loads.
Journal of Sound and Vibration., 180(4) :637644, 1995.
[49] K. Henchi, M. Fafard, G. Dhatt, and M. Talbot. Dynamic behaviour of multi-span beams
under moving loads. Journal of Sound and Vibration, 199(1) :3350, 1997.
[50] H.M. Hilber, T.J.R. Hughes, and R.L. Taylor. Improved numerical dissipation for time inte-
gration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics,
5 :283292, 1977.
[51] T.J. Hughes. The Finite Element Method - Linear Static and Dynamic Finite Element
Analysis. Dover Publications, 2000.
[52] T.J.R. Hughes and G.M. Hulbert. Space-time element methods for elastodynamics : for-
mulation and error analysis. Computer Methods in Applied Mechanics and Engineering,
66 :339363, 1987.
[53] G.M. Hulbert. Time finite element methods structural dynamics. International Journal for
Numerical Methods in Engineering, 33 :307331, 1992.
[54] H.H. Hung and Y.B. Yang. Elastic waves in visco-elastic half-space generated by various
vehicle loads. Soil Dynamics and Earthquake Engineering., 21 :117, 2001.
[55] IWRN-5. Fifth International Workshop on Railway an Tracked Transit System. Journal of
Sound and Vibration (Special issue), 193(1), 1996.
[56] K.L. Johnson. Contact mechanics. Cambridge University Press, 1992.
[57] C.J.C. Jones, X. Sheng, and M. Petyt. Simulation of ground vibration from a moving har-
monic load on a railway track. Journal of Sound and Vibration, 231(3) :739751, 2000.
[58] D.V. Jones, D. Le Houedec, A.T. Peplow, and M. Petyt. Ground vibration in the vicinity of
a moving harmonic rectangular load on a half-space. Eur. J. Mech. A/Solids, 17(1) :153166,
1998.
[59] D.V. Jones and M. Petyt. Ground vibration in the vicinity of a strip load : an elastic layer
on an elastic half-space. Journal of Sound and Vibration, 61(1) :118, 1993.
172 BIBLIOGRAPHIE
[60] D.V. Jones and M. Petyt. Ground vibration in the vicinity of a vicinity of a rectangular load
on a half-space. Journal of Sound and Vibration, 161(1) :141159, 1993.
[61] P.H. Kirkegaard, S.R.K. Nielsen, S. Krenk, and L. Kellezi. Radiation of Air-Borne Noise in
Non-Homogeneous Wind and Temperature Fields using FEM Analysis. In Fryba & Naprstek,
editor, Structural Dynamics - EURODYN 99, Rotterdam, 1999. Balkena.
[62] K. Knothe and Y. Wu. Receptance behaviour of railway track and subgrade. Archive of
Applied Mechanics, 68 :457470, 1998.
[63] A.W.M Kok. Lumped impulses, discrete displacements and moving load analysis. HERON,
42(1), 1997.
[64] A.W.M. Kok. Finite element models for the steady state analysis of moving loads. HERON,
45(1) :5361, 2000.
[65] S. Krenk, L. Kellezi, S.R.K. Nielsen, and P.H. Kirkegaard. Finite elements and transmitting
boundary conditions for moving loads. In Fryba & Naprstek, editor, Structural Dynamics -
EURODYN 99, Rotterdam, 1999. Balkena.
[66] V. Krylov. Generation of ground vibration by superfast trains. Applied Acoustics, 44(149
164), 1995.
[67] V. Krylov. Vibration impact of high-speed trains. I. Effect of track dynamics. The Journal
of the Acoustical Society of America, 100 :31213134, 1996.
[68] V. Krylov. Vibration impact of high-speed trains. I. Effect of track dynamics. The Journal
of the Acoustical Society of America, 101 :3810, 1997. Erratum.
[69] V. Krylov, A.R. Dawson, M.E. Heelis, and A.C. Collop. Rail movement and ground waves
caused by high-speed trains approaching track-soil critical velocities. In Pro. Instn. Mech.
Engrs., volume 214, pages 107116, 2000.
[70] T. Krzyzynski, T. Popp, and H. Kruse. The dynamics of raiway tracks based on the theory
of periodic structures. ZAMM - Z. Angew. Math. Mech., 78, 1998. S1.
[71] P. Le Tallec. Introduction a la dynamique des structures. Edition Ellipses, 2000. (Cours de
lEcole Polytechnique).
[72] G. Lefeuvre-Mesgouez and D. Le Houedec. Propagation dondes dans un massif soumis a
des charges roulantes se deplacant a vitesse elevee. In 14e Congres Francais de Mecanique,
Toulouse, 1999.
[73] G. Leufeuve-Mesgouez. Propagation d ondes dans un massif soumis a des charges se
deplacant a vitesse constante. PhD thesis, Ecole Centrale de Nantes, 1999.
[74] X.D. Li and N.-E. Wiberg. Implementation and adaptivity of a space-time finite element
method for structural dynamics. Computer Methods in Applied Mechanics and Engineering,
156 :211229, 1998.
[75] M. Lieb and B. Sudret. A fast algorithm for soil dynamics calculations by wavelet decompo-
sition. Archive of Applied Mechanics., 68 :147157, 1998.
[76] G. Lombaert, G. Degrande, and D. Clouteau. Numerical modelling of free field traffic-induced
vibration. Soil Dynamics and Earthquake Engineering., 19(473488), 2000.
[77] M. Lucchesi, C. Padovani, and G. Pasquinelli. On the numerical solution of equilibrium
problems for elastic solids with bounded tensile strength. Comput. Methods Appl. Mech.
Engrg., 127 :3756, 1995.
[78] F. Maddalena and M. Ferrari. Viscoelasticity of granular materials. Mechanics of Materials,
20(241250), 1995.
[79] S. Melin. Wave propagation in granular assemblie. Physical Review, 49(3), 1994.
BIBLIOGRAPHIE 173
[80] A. Metrikine and H. Dieterman. Instability of vibration of a mass moving uniformly along
an axial compressed beam on a viscoelastic foundation. Journal of Sound and Vibration,
201(5) :567576, 1997.
[81] A. Metrikine and K. Popp. Vibration of a periodically supported beam on an elastic half-
space. Eur. J. Mech., A/Solids, 18 :679701, 1999.
[82] A. Metrikine and A.C.M. Vrouwenvelder. Surface ground vibration due to a moving train in
a tunel : two-dimensional model. Journal of Sound and Vibration, 234(1) :4366, 2000.
[83] A.V. Metrikine and H.A. Dieterman. Lateral vibration of an axially compressed beam on an
elastic half-space due to a moving lateral load. Eur. J. Mech. A/Solids, 18 :147158, 1999.
[84] H.-B Muhlhaus and F.. Oka. Dispersion and wave propagation in discrete and continuous
models for granular materials. International Journal in Solids & Structures, 33 :28412858,
1996.
[85] V.H. Nguyen, D. Duhamel, and B. Nedjar. Calcul de structures soumises a une charge
mobile par la methode des elements finis. XVeme Congres Francais de Mecanique, Nancy,
37 Septembre, 2001.
[86] V.H. Nguyen, D. Duhamel, and B. Nedjar. A no-tension model for granular material. In 2nd
European Conference on Computational Mechanics - ECCM2001, June 2001.
[87] M. Olsson. On the fundamental moving load problem. Journal of Sound and Vibration,
145(2) :299307, 1991.
[88] X. Oviedo. Etude du comportement du ballast par un modele micromecanique. PhD thesis,
Ecole Nationale des Ponts et Chaussees, May 2001.
[89] X. Oviedo, K. Sab, and P.-E. Gautier. Etude du comportement du ballast ferroviaire par
un modele micromecanique. Colloque Physique et Mecanique des Materiaux Granulaires.
Champs-sur-Marne, 5-7 Septembre, 2000. Actes des journees scientifiques du LCPC, tome
2, pages 339344, 2000.
[90] G. Pan, H. Okada, and S.N. Atluri. Nonlinear transient dynamic analysis of sol-pavement
interaction under moving load : a couple BEM-FEM approach. Engineering Analysis with
Boundary Element., 14(99112), 1994.
[91] T. Panzeca and G. Del Piero. Constitutive equations for no-tension materials. Meccanica,
23 :8893, 1988.
[92] B. Picoux and D. Le Houedec. Propagation dondes a la surface dun massif sousmis a un
effort lateral et vertical du a un traffic ferroviaire. In XVeme Congres Francais de Mecanique,
Nancy, 2001.
[93] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C -
The Art of Scientific Computing. Cambridge University Press, second edition, 1992.
[94] V. Profillidis. La voie ferree et sa fondation - modelisation mathematique. PhD thesis, Ecole
Nationale des Ponts et Chaussees, 1983.
[95] K.M. Rassmusen, S.R.K. Nielsen, and P.H. Kirkegaard. Boundary element method solution
in the time domain for a moving time-dependent force. Computers & Structures, 79 :691701,
2001.
[96] C. Renard and P. Orsero. Influence de limperfection de la voie sur le comportement dy-
namique du pont au passage dun train. 14ieme Congres Fracais de Mecanique, Toulouse,
1999.
[97] J. R. Rieker, Y.-H. Lin, and M.W. Trethewey. Discretisation consideration in moving load
finite element beam models. Finite Element in Analysis and Design, 21 :129144, 1996.
[98] M.H. Sadd, J. Gao, and A. Shukla. Numerical analysis of wave propagation through assem-
blies of elliptical particles. Computers and Geotechnics, 20(3/4) :323343, 1997.
174 BIBLIOGRAPHIE
[99] R. Sauvage and G. Larible. La modelisation par elements finis des couches dassise de la
voie ferree. Revue Generale des Chemins de Fer., September 1982.
[100] X. Sheng, C.J.C. Jones, and M. Petyt. Ground vibration generated by a harmonic load acting
on a railway track. Journal of Sound and Vibration, 225(1) :328, 1999.
[101] X. Sheng, C.J.C. Jones, and M. Petyt. Ground vibration generated by a load moving along
a railway track. Journal of Sound and Vibration, 228(1) :129156, 1999.
[102] A. Shulka, M.H. Sadd, and H. Mei. Experimental and computational modelling of wave
propagation in granular materials. Experimental Mechanics, 30(3), 1991.
[103] A.S.J. Suiker, C.S. Chang, R. de Borst, and C. Esveld. Surface waves in a stratified half space
with enhanced continuum properties. Part 1 : Formulation of the boundary value problem.
European Journal of Mechanics. A/Solids, 18(749768), 1999.
[104] A.S.J. Suiker, C.S. Chang, R. de Borst, and C. Esveld. Surface waves in a stratified half
space with enhanced continuum properties. Part 2 : Analysis of the wave characteristics in
regard to high-speed railways track. European Journal of Mechanics. A/Solids, 18(769784),
1999.
[105] A.S.J. Suiker, A.V. Metrikine, and R. De Borst. Dynamic behaviour of a layer of discrete
particles. Part 2 : Response to a uniformly moving harmonical vibrating load. Journal of
Sound and Vibration., 240(1) :1939, 2001.
[106] C. Wassgren. Vibration of Granular Materials. PhD thesis, Institut Califonia of Technology,
1997.
[107] J.P. Wolf and C. Song. Finite-Element Modelling of Unbounded Media. Wiley, NewYork,
1996.
[108] A.R.M. Wolfert, H.A. Dietermann, and A.V. Metrikine. Passing through the elastic wave
barrier by a load moving along a waveguide. Journal of Sound and Vibration, 203(4) :597
606, 1997.
[109] J.-J. Wu, A.R. Whittaker, and M.P. Cartmell. The use of finite element techniques for
calculating the dynamic response of structures to moving loads. Computers anf Structures,
78 :789799, 2000.
[110] Y.-B. Yang and H.-H. Hung. A 2.5D finite/infinite element approach for modelling visco-
elastic bodies subjected to moving loads. International Journal for Numerical Methods in
Engineering., 51(13171336), 2001.
[111] H.R. Yerli, B. Temel, and E. Kiral. Multi-wave transient and harmonic infinite elements for
two-dimensional unbounded domain problems. Computers and Geotechnics, 24(185206),
1999.
[112] Y.-G. Zhang and J. Ballman. Two techniques for the absorption of elastic waves using an
artificial transition layer. Wave Motion, 25(1533), 1997.
[113] Y. Zhu, A. Shukla, and M.H. Sadd. The effect of microstructural fabric on dynamic load trans-
fer in two dimensional assemblies of elliptical particles. J. Mech. Phys. Solids, 44(8) :1283
1303, 1996.
[114] H.S. Zibdeh and R Rackwitz. Moving loads on beams with general boundary conditions.
Journal of Sound and Vibration., 195(1) :85102, 1996.
[115] O.C. Zienkiewicz and R.L. Taylor. The Finite Element Method. Volume 2 : Solid and Fluid
Mechanics Dynamics and Non-linearity. McGraw-Hill, 1991.
Resume
Comportement dynamique de structures non-lin eares soumises a des charges mobiles
Ce travail a pour but letude de la dynamique de structures non-lineaires soumises a des passages de
vehicules. Le cas dun demi-espace multi-couche visco-elastique soumis a une charge mobile est premierement
resolu par une methode semi-analytique. Il permet de comprendre les phenomenes physiques et peut servir
a valider des approches plus complexes elaborees dans la suite. Plusieurs etudes sur ce probleme ont ete
realisees mais il en existe encore peu qui tiennent compte de la non-linearite des materiaux. Ces non-
linearites, qui peuvent etre tres forte dans certains materiaux (par exemple le ballast), causent des surcharges
dynamiques plus importantes. On considere ici deux types de non-linearite : lelasticite non-lineaire et
lunilateralite. La structure est etudiee avec des modeles tres simples en unidimensionnel ou plus complexes
en tridimensionnel. La methode des elements finis est choisie pour les calculs. Le cas 1D consiste en un
probleme uniaxial dune barre et un probleme de flexion dune poutre posee sur un systeme de ressorts-
amortisseurs. Le cas 3D utilise un modele de materiau ballast qui est elastique (non-lineaire) en compression
et unilateral dans les directions principales. Comme on sinteresse au regime permanent, on ecrit les equations
des elements finis et on les resout dans le repere mobile en utilisant un changement de variable. Dans le cas
ou la vitesse de la charge devient supersonique, la methode des elements finis dans le repere mobile nest
plus valable et on doit utiliser une autre approche dans laquelle le changement de variable est realise a partir
de lequation discretisee. On montre aussi pour le modele uniaxial, dans le cas ou la vitesse de la charge
se situe entre la vitesse de londe de compression et la vitesse de londe de tension, quil peut exister des
ondes de chocs dans le regime permanent. Les resultats obtenus sont presentes et sont compares avec ceux
des cas lineaires. Ils montrent les influences importantes de la non-linearite sur la reponse dynamique des
structures.
Mots cles : charge mobile, comportement unilateral, dynamique nonlineaire, elements finis non-lineaires,
regime permanent.
Abstract
Dynamical behaviour of nonlinear structures under moving loads
The aim of this work is to study the nonlinear dynamic of structures under passages of vehicles. The
problem of a multi-layer viscoelastic halfspace under a moving load is firstly solved by a semi-analytical
method. It allows to understand the physical phenomena and may be used to valid the more complex
approaches established later. Many studies on this problem have been realised but still few of its take
into account the non-linearity of materials. These non-linearities, which may be very important in some
materials (e.g. the ballast), cause the more important dynamical overloadings. We consider here two kinds
of nonlinearity : the nonlinear elastic behaviour and the unilaterality. The structure is studied with simple
models in one dimension or more complex in three dimensions. The finite element method is chosen in
the analyses. The 1D case involves the problem of an uniaxial bar and the problem of a flexion beam
posed on viscoelastic strings system. The 3D case uses a ballast material model which is (nonlinear) elastic
in compression and is unilateral in principal directions. As we are interested in the stationnary solution,
we establish and resolve the finite element equations in the moving coordinates by using of a variable
transformation. When the load velocity becomes supersonic, the finite element method is no longer valid
and we have to use another approach in which the variable transformation is based on discretizated equations.
We show also for the uniaxial model, when the load velocity is between the pression wave velocity and the
tension ones, that the shock waves may appear in the stationnary state. The obtained results are presented
and compared with the linear analysiss ones. These results show the important influences of the nonlinearity
on the dynamical response of structures.
Keywords : moving load, unilateral constitutive, nonlinear dynamic, nonlinear finite element, stationnary
state.