0% ont trouvé ce document utile (0 vote)

1 vues46 pages

Chap 1

Ce document présente une étude bibliographique sur la mécanique de la rupture, abordant les concepts de base, les critères de propagation, et les méthodes de modélisation par éléments finis. Il décrit également les différents types de rupture des matériaux, leur classification, et les mécanismes sous-jacents. Enfin, il explore les approches théoriques et numériques pour la propagation de fissures dans les structures.

Transféré par

bikoyelieemmanuel97

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

1 vues46 pages

Chap 1

Transféré par

bikoyelieemmanuel97

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

Chapitre 1

Étude bibliographique

Sommaire
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Mécanique de la rupture . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.1 Concepts de base et définitions . . . . . . . . . . . . . . . . . . . . 10
1.2.2 Critères de propagation . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.3 Critères de direction de propagation . . . . . . . . . . . . . . . . . 19
1.2.4 Vitesse de propagation . . . . . . . . . . . . . . . . . . . . . . . . 21
1.3 Modélisation éléments finis d’une structure mince . . . . . . . . . . . . 23
1.3.1 Discrétisation par la méthode des éléments finis . . . . . . . . . . . 23
1.3.2 Cadre d’utilisation des coques . . . . . . . . . . . . . . . . . . . . 26
1.3.3 Théorie de Kirchhoff-Love . . . . . . . . . . . . . . . . . . . . . . 31
1.3.4 Théorie de Mindlin-Reissner . . . . . . . . . . . . . . . . . . . . . 32
1.3.5 Verrouillage numérique . . . . . . . . . . . . . . . . . . . . . . . . 34
1.3.6 Quelques éléments finis de coques classiques . . . . . . . . . . . . 37
1.4 Modélisation numérique de la propagation de fissure . . . . . . . . . . 38
1.4.1 Méthodes non basées sur les éléments finis . . . . . . . . . . . . . 38
1.4.2 Méthodes basées sur les éléments finis classiques . . . . . . . . . . 40
1.4.3 Méthodes basées sur la partition de l’unité . . . . . . . . . . . . . . 45
1.4.4 Méthodes diffuses . . . . . . . . . . . . . . . . . . . . . . . . . . 48
1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Cette thèse est accessible à l'adresse : [Link]

1.1 Introduction
Après une brève introduction sur les différents types de rupture qui peuvent inter-
venir en fonction du matériau et du chargement considéré, les principes généraux de la
mécanique de la rupture qui permettent d’accéder aux différents critères de propagation
employés dans ces travaux sont exposés. Ensuite sont abordés les différents aspects de la
modélisation par éléments finis d’une structure mince. Enfin, une dernière partie présente
les méthodes numériques qui existent actuellement pour représenter puis faire propager
une fissure. On décrira notamment la méthode des éléments finis étendus.

Qu’est-ce que la rupture ?

Les mécanismes de rupture observés sont très variés et fortement dépendant du maté-
riau et de l’essai comme on peut le voir figure 1.1. Par exemple lors d’un essai de traction
sur acier, on va d’abord observer un phénomène dit de striction qui se caractérise par
un rétrécissement de la section de l’éprouvette. La rupture intervient ensuite par une sé-
paration en deux de l’éprouvette dans le sens transversal. Pour un essai de compression
d’un composite dans la direction des fibres, la rupture a pour origine une décohésion de
celles-ci. Pour un béton en compression, la rupture se traduit par l’apparition d’un réseau
de fissures longitudinales alors que dans un essai de traction elle consiste essentiellement
en une fissure unique transversale.

F IGURE 1.1: État final de différentes éprouvettes après un essai uniaxial (tiré de
[MAR 12]).

Une classification du comportement et de la rupture des matériaux peut être proposée

à partir des différences observées dans des essais uniaxiaux. Lorsqu’un matériau comme
le béton voit sa contrainte rapidement chuter vers 0 après la phase élastique, il est qualifié
de fragile. Alors qu’un matériau comme l’acier à température ambiante peut supporter
des déformations importantes avant de rompre, il est dit ductile. Entre les deux, on

Cette thèse est accessible à l'adresse : [Link]

trouve des matériaux dits quasi-fragiles pour lesquels les déformations plastiques restent
confinées à une zone restreinte proche de la singularité.

F IGURE 1.2: Réponses schématiques de différents matériaux dans un essai de

traction uniaxiale.

Cependant cette distinction fragile-ductile ne permet ni de caractériser un matériau

donné, ni de classer l’ensemble des matériaux car leurs réponses peuvent fortement
dépendre de deux paramètres : la température et la vitesse de déformation. Ainsi, un
métal qui est ductile à température ambiante peut devenir fragile à basse température.
Cette transition fragile-ductile, longtemps méconnue ou mésestimée est d’ailleurs à la
base de plusieurs accidents dont les célèbres Liberty ships (voir F IG .1.3). Ces cargos
construits durant la Seconde Guerre mondiale qui naviguaient dans l’Atlantique Nord
étaient exposés à des températures inférieures à la température de transition, faisant
passer le matériau de la coque d’un état ductile à un état fragile à cause d’une chute
sensible de la ténacité. Au total, plus de 200 navires sur les 5000 construits ont subit des
dégâts importants.

Les différents types de ruptures peuvent être classés selon deux méthodes distinctes :
la première prend en compte à leur cinétique (F IG .1.4) et la seconde s’appuie sur les
observations microstructures (F IG .1.5).

Les ruptures ductiles sont toujours caractérisées par la présence d’une importante zone
déformée plastiquement. Elles font apparaître une surface matte en relief accompagnée
d’une zone de striction si la sollicitation est de traction. On ne distingue pas de zone
d’amorçage précise ni de direction de propagation nette. L’analyse microscopique du fa-
ciès révèle de nombreuses cupules et éventuellement des inclusions. Ces ruptures sont
transgranulaires (voir F IG .1.6).

Les ruptures fragiles sont exemptes de déformations plastiques. Leur surface est plane
et brillante. On ne distingue ni zone d’amorçage ni direction de propagation. Leur origine

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.3: Photographie d’un Liberty ship dont la coque est rompue en deux.

F IGURE 1.4: Classement des ruptures selon leur cinétique.

F IGURE 1.5: Classement des ruptures selon leur observation.

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.6: Faciès d’une rupture ductile avec présence de cupules (à gauche) et
d’inclusions (à droite). (tiré de [LOU 14]).

peut trouver sa source dans plusieurs mécanismes physiques. Le mécanisme de clivage

est lui aussi transgranulaire. À l’échelle microscopique, on note la présence de plans et de
rivières de clivage comme on peut le voir figure 1.7.

F IGURE 1.7: Faciès d’une rupture fragile par clivage (tiré de [LOU 14]).

Les mécanismes de décohésion présentent un aspect globalement plan et ne font pas

apparaître de zone d’amorçage et de direction de propagation. À l’échelle microscopique,
elles font apparaître des polyèdres juxtaposés correspondant à la surface des grains défi-
nissant le faciès de rupture. Cela est dû au fait que ces ruptures sont intergranulaires. Des
phases fragilisantes peuvent apparaître au niveau de la surface de rupture. Bien souvent,
des fissures secondaires sont observées (voir F IG .1.8).

F IGURE 1.8: Rupture fragile par décohésion (tiré de [LOU 14]).

Cette thèse est accessible à l'adresse : [Link]

1.2 Mécanique de la rupture

Cette partie présente les notions théoriques de la mécanique de la rupture. Après un
bref historique, les hypothèses de bases sont rappelées puis les différentes approches pour
traiter la propagation de fissure sont exposées.

1.2.1 Concepts de base et définitions

La mécanique de la rupture a pour objet d’étudier et de prédire l’amorçage et la pro-
pagation de fissures. Seule, elle ne permet pas de modéliser le scénario entier de ruine de
la structure avec notamment la phase de localisation. En revanche, à partir d’un défaut
existant, la mécanique de la rupture permet d’obtenir des résultats satisfaisant avec une
modélisation grossière.
Historiquement, cette branche de la mécanique s’est fortement développée au début
du XX e siècle alors que les industriels tentaient de prévoir si les défauts présents dans les
structures présentaient des risques de propagation de fissures.

La rupture est caractérisée par la séparation irréversible d’un milieu continu en deux
parties de part et d’autre d’une interface A (F IG .1.9).

F IGURE 1.9: Séparation d’un milieu de part et d’autre d’une Interface A.

Cette séparation est appelée fissure et modifie les champs de déplacements, de défor-
mations et de contraintes. La discontinuité du champ de déplacement se caractérise par
−
un saut de déplacement [ui ] = u+ i − ui . D’un point de vue cinématique, on peut définir
trois modes de rupture comme présenté sur la figure 1.10 :
– le mode I ou mode d’ouverture : le déplacement est perpendiculaire au plan de la
fissure,
– le mode II ou mode de cisaillement plan : le déplacement est parallèle au plan de la
fissure et normal au front,
– le mode III ou mode de cisaillement antiplan : le déplacement est parallèle au plan
de la fissure et au front.

Cette thèse est accessible à l'adresse : [Link]

Bien souvent, la fissure se propage suivant une combinaison de ces modes. On parle alors
de mixte.

F IGURE 1.10: Modes de rupture.

Inglis [ING 13] est le premier à réaliser des études pour évaluer les concentrations
de contraintes dans le voisinage d’une pointe de fissure. Il a notamment résolu analyti-
quement le cas élastique d’un trou elliptique dans une plaque soumise à un chargement
uniforme. En faisant tendre le rayon de courbure vers zéro, ses résultats lui ont permis de
conclure à une contrainte infinie en pointe d’une fissure très fine. Ces travaux sont à la
base de l’analyse asymptotique.

[Link] Approche énergétique et analyse asymptotique

C’est cette divergence du champ contrainte en pointe du fissure qui a motivé Griffith
[GRI 21] à entreprendre une approche énergétique de la mécanique de la rupture. Il fait
notamment le lien entre la perte d’énergie du système et l’avancée de la fissure en élasticité
linéaire.
L’énergie ainsi dépensée pour permettre l’accroissement d’une fissure d’une aire dA
est appelée taux de restitution de l’énergie et noté G en l’honneur de Griffith. Le compor-
tement de la fissures se caractérise donc par le transfert de l’énergie potentielle E p de la
structure en énergie de décohésion caractérisée par γ au voisinage de la pointe.

dE p
G=− = 2 γ. (1.1)
dA
L’analyse locale en élasticité linéaire permet de montrer l’existence d’une solution
asymptotique singulière pour les champs mécaniques. Irwin a alors proposé de définir des
facteurs quantifiant l’intensité de la singularité des contraintes [IRW 56]. Ainsi pour un

Cette thèse est accessible à l'adresse : [Link]

chargement statique (ou quasi-statique) ces facteurs sont utilisés pour déterminer l’inten-
sité de la singularité aussi bien en terme de contraintes que de déplacements. Ces facteurs
sont appelés facteurs d’intensité des contraintes et sont définis pour chaque mode élé-
mentaire de rupture :

√
r
µ 2π
KI = lim 2πr σ22 (θ = 0) = lim Ju2 (θ = π)K, (1.2)
r→0 r→0 k + 1 r
√
r
µ 2π
KII = lim 2πr σ12 (θ = 0) = lim Ju1 (θ = π)K, (1.3)
r→0 r→0 k + 1 r
√
r
µ 2π
KIII = lim 2πr σ23 (θ = 0) = lim Ju3 (θ = π)K, (1.4)
r→0 r→0 4 r

où Jui (θ = π)K = ui (θ = π) − ui (θ = −π) sont les sauts de déplacements au passage de la

discontinuité et k la constante de Kolosov, telle que :
(
3 − 4ν en déformations planes,
k = 3−ν (1.5)
1−ν en contraintes planes.

F IGURE 1.11: Repère local en pointe de fissure.

On peut également introduire ces facteurs dans l’expression des solutions analy-
tiques en pointe de fissure. En élasticité plane, l’approche de Westergaard reprise par
Bui [BUI 78] permet d’obtenir les déplacements et les contraintes au voisinage du front à
l’aide des fonctions d’Airy :
r
1 r θ θ
u1 (r, θ) = KI cos (k − cosθ) + KII sin (k + cos θ + 2) , (1.6)
2µ 2π 2 2
r
1 r θ θ
u2 (r, θ) = KI sin (k − cosθ) − KII cos (k + cos θ − 2) , (1.7)
2µ 2π 2 2
r
2 r θ
u3 (r, θ) = KIII sin . (1.8)
µ 2π 2

Cette thèse est accessible à l'adresse : [Link]

L’expression des contraintes exprimées dans le repère cylindrique centré en pointe de

fissure donne :

KI θ 3θ KII θ 3θ
σrr = √ 5 cos − cos − √ 5 sin − 3 sin , (1.9)
4 2πr 2 2 4 2πr 2 2

K θ 3θ KII θ 3θ
σθθ = √I 3 cos + cos − √ 3 sin + 3 sin , (1.10)
4 2πr 2 2 4 2πr 2 2

KI θ 3θ KII θ 3θ
σrθ = √ sin + sin + √ cos + 3 cos , (1.11)
4 2πr 2 2 4 2πr 2 2

K θ
σrz = √ III sin , (1.12)
2πr 2

KIII θ
σθz = √ cos , (1.13)
2πr 2
(
ν(σrr + σθθ ) en déformation plane,
σzz = (1.14)
0 en contrainte plane,

Irwin a également permis de directement relier les facteurs d’intensité des contraintes
au taux de restitution de l’énergie :
2
KI2 + KII2 KIII
G= + , (1.15)
E∗ 2µ
(
E
E 2 en déformations planes,
avec µ= et E ∗ = 1−ν
2(1 + ν) E en contraintes planes.
Limités jusqu’alors aux ruptures des matériaux fragiles, Irwin [Irw60] a ensuite
étendu le concept de Griffith aux matériaux quasi fragiles en utilisant la plasticité pour
majorer la valeur du champ de contrainte en pointe de fissure.

En 1968, Rice [RIC 68b] a introduit le concept de l’intégrale J permettant d’évaluer le

taux de restitution de l’énergie G. L’intégrale J utilise un contour fermé Γ qui entoure la
pointe de fissure, comme illustré figure 1.12. Cette intégrale est un invariant indépendant
du contour Γ choisi, ce qui est un résultat important pour les méthodes numériques. Cela
permet d’évaluer G avec des informations situées loin de la pointe de fissure qui perturbe
sensiblement la qualité des champs calculés dans son voisinage.

Le passage à la dynamique, avec la prise en compte des effets d’inertie tant sur le
mouvement de la structure que sur celui de la fissure est réalisé par Craggs [CRA 60] et
par Freund [FRE 72]. Ces travaux ont permis de généraliser les concepts de la mécanique
de la rupture à des phénomènes transitoires. On peut parler de rupture dynamique lorsque
le chargement est rapide comme c’est le cas lors d’un crash, ou bien lorsque la vitesse
de propagation de la fissure est non négligeable par rapport aux vitesses des ondes dans
le matériau. Lors d’une propagation dynamique, il devient nécessaire de différencier le
facteur d’intensité des contraintes et le facteur d’intensité des déplacements.

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.12: Contour pour l’intégrale J.

1.2.2 Critères de propagation

[Link] Mécanique linéaire de la rupture
Pour la propagation d’une fissure donnée, il existe une valeur critique du facteur d’in-
tensité des contraintes notée Kc qui est directement reliée à la ténacité du matériau. La
ténacité√est la capacité d’un matériau à résister à la propagation d’une fissure et se mesure
en MPa m. On peut la définir définir comme étant la quantité d’énergie qu’un matériau
peut absorber avant de rompre. C’est une propriété matériau indépendante de la géométrie
de l’éprouvette qui se détermine expérimentalement. Ainsi, pour un essai sollicitant une
éprouvette en mode I, il y a initiation de la fissure lorsque KI = KIc .
De la même manière, on peut définir un taux de restitution de l’énergie critique, noté
Gc = 2γ, qui s’apparente à un seuil pour l’initiation de la propagation de fissure. Le
critère de propagation de Griffith est donc formulé de la manière suivante :

– G < Gc absence de propagation,

– G = Gc initiation ou propagation stable,
– G > Gc propagation instable.

Ces deux approches en taux de restitution et en intensité des contraintes sont équi-
valent. On peut passer de l’un à l’autre par la formule suivante en mode I :
(KIc )2
Gc = . (1.16)
E∗
Dans le cas des modes mixtes, la relation devient :
(KIc )2 + (KIIc )2 (KIIIc )2
Gc = + , (1.17)
E∗ 2µ
(
E
E 2 en déformations planes,
avec µ= et E ∗ = 1−ν
2(1 + ν) E en contraintes planes.

Cette thèse est accessible à l'adresse : [Link]

[Link] Mécanique élasto-plastique de la rupture

Pour les matériaux au comportement fortement non linéaire tels que les matériaux
ductiles, la mécanique linéaire élastique de la rupture ne suffit plus et il faut alors se
placer dans le cadre de la mécanique élasto-plastique de la rupture. Différents modèles
permettent de prévoir la propagation d’une fissure

Crack Tip Opening Displacement (CTOD) Wells [WEL 61] puis Cotterel [COT 62]
ont introduit la notion de CTOD (crack-tip opening displacement) qui établit un critère de
propagation basé sur l’ouverture en pointe de fissure pour des matériaux élasto-plastiques.
Lorsque la contrainte en pointe de fissure atteint une valeur critique, la fissure commence
à propager. La CTOD notée δ est décrite comme cette séparation des lèvres de la fissure.
Lorsqu’elle atteint une valeur critique notée δc , la fissure propage.

Intégrale-J Rice [RIC 68a, RIC 68b] a lui adapté le concept d’intégrale J afin d’évaluer
le taux de restitution de l’énergie pour les matériaux non linéaires. Contrairement au
cas linéaire élastique, la séparation des lèvres de la fissure n’est plus la seule cause de
dissipation d’énergie. Une partie de l’énergie est en effet dépensée dans le phénomène de
plasticité. L’idée de Rice est de considérer le matériau non pas comme élasto-plastique
mais élastique non linéaire. Comme on peut le voir figure 1.13, les deux types de com-
portement sont identiques tant qu’il n’y a pas de décharge. Le matériau élasto-plastique
laisse en effet une déformation plastique qui est non réversible si le chargement est retiré.

F IGURE 1.13: Courbe de traction d’un matériau élastique non linéaire et d’un
matériau élasto plastique.

Ainsi, l’intégrale J proposée par Rice reste valide et indépendant du contour Γ tant
qu’il n’y a pas de décharge. Lorsque la fissure a propagé cette hypothèse de non décharge

Cette thèse est accessible à l'adresse : [Link]

n’est plus vérifiée en arrière de la pointe. En avant, cette hypothèse reste raisonnable. Il
en résulte que l’indépendance au contour n’est plus garantie.

Approche locale en pointe de fissure D’autre critères basés sur l’analyse des champs
mécaniques locaux existent. On peut citer les travaux de Ritchie, Knott et Rice sur la
rupture par clivage [RIT 73], basés sur l’observation des grandeurs mécaniques à une
distance rc en avant de la pointe de fissure. Il y a alors propagation si la contrainte locale
d’ouverture en ce point amont est supérieure à une contrainte critique. Il a été montré que
ce critère RKR n’était pas complètement satisfaisant car très dépendant de la distance rc .
De plus, les contraintes locales oscillent fortement en dynamique rendant le critère plus
"aléatoire".
Jirásek [JIR 98] propose de moyenner les grandeurs en pointe de fissure. L’idée n’est
pas de délocaliser le tenseur des contraintes comme dans l’approche non-locale de type
Pijaudier-Cabot et Bazant [PIJ 87] pour l’endommagement, mais de s’affranchir de pos-
sibles erreurs numériques locales en intégrant les contraintes sur un petite zone et en les
pondérant d’une fonction poids qui dépend de la distance à la pointe de fissure. Dans les
travaux de Wells [WEL 01] et Remmers [REM 03] cette zone est un disque en pointe de
fissure de rayon égal à trois fois la taille d’un élément. La fonction de pondération choisie
est de la forme d’une Gaussienne w et est définie par l’équation 1.18 (voir F IG .1.14) :

r 2
w(r) = exp − . (1.18)
R

Le tenseur des contraintes moyenné σ̃i j se calcule donc avec l’équation 1.19 :

−( Rr )2 dS
R
D σi j (M)e M
σ̃i j = R −( r )2 , (1.19)
D e R dSM

Ainsi, plus un point d’intégration (lieu de calcul de la contrainte) est proche de la

pointe de fissure, plus le poids qui lui est associé est important. La distance R sur laquelle
le tenseur des contraintes est moyenné introduit une notion de distance caractéristique
comme le critère RKR. Elle doit être choisie ni trop grande pour suffisamment caractériser
l’état en pointe de fissure, ni trop faible pour définir une zone de dimension raisonnable
par rapport à la taille d’un élément. En général on prend R est pris égal à un nombre entier
d’éléments.
Le critère est ensuite formulé en considérant qu’il y a propagation lorsque la contrainte
principale maximale associée à σ˜i j est supérieure à une valeur critique notée σIc . L’intérêt
de cette méthode est qu’elle est utilisable à la fois pour un comportement linéaire ou non
linéaire du matériau, mais elle nécessite d’identifier σIc . Des travaux récents [PRA 07,
MEN 07, HAB 12b, SAI 14] sont basés sur l’utilisation de ce critère.

Fracture Process Zone Enfin, le caractère non linéaire du matériau proche de la pointe
de fissure peut être modélisé en introduisant la notion de Fracture Process Zone (FPZ).

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.14: Fonction poids w(r) en fonction du rapport r/R.

Cette zone peut être définie comme le lieu où les défauts croissent et coalescent pour
former la fissure dans le cas d’un matériau ductile. Plutôt que de décrire ce phénomène
avec la mécanique de l’endommagement, on peut ajouter au modèle une cohésion entre
les lèvres de fissures à l’aide des modèles de zones cohésives introduits par Dugdale
[DUG 60] et Barenblatt [BAR 62].
Cette méthode fut d’abord appliquée avec les éléments finis pour réaliser des
calculs de propagation de fissures dans le béton [HIL 76], puis de la propagation de
fissures dans les matériaux ductiles [NEE 87]. Cette approche est intermédiaire entre les
modèles basés sur la mécanique de la rupture utilisant des critères d’amorçage (taux de
restitution d’énergie ou critères en contrainte) et ceux utilisant la mécanique continue
de l’endommagement. Lorsque la localisation géométrique des fissures potentielles est
connue, elle fournit un cadre sûr et robuste pour prédire la cinétique de propagation.

Dans la process zone, une relation dite de traction séparation existe entre les
contraintes cohésives et le saut de déplacement entre les lèvres de la fissure. Physique-
ment, l’aire sous la courbe (F IG .1.15) représente le travail dissipé par l’ouverture des
lèvres identifié comme le taux de restitution d’énergie. Une fissure est alors définie par
trois zones, la première correspond à une fissure libre de contrainte, la seconde, appelée
zone cohésive, où les forces de cohésion sont non nulles et enfin une zone saine (ou zone
de fissuration potentielle) où la contrainte d’amorçage n’a pas été atteinte (voir F IG .1.16).

Le critère de propagation se fait alors en énergie : il y a propagation lorsque l’énergie

apportée est supérieure à celle décrite par la loi de traction-séparation (voir F IG .1.15)
injectée dans la zone cohésive.

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.15: Loi de traction séparation.

F IGURE 1.16: Schématisation des zones cohésives.

Cette thèse est accessible à l'adresse : [Link]

1.2.3 Critères de direction de propagation

Pour simuler la propagation de fissure en se basant sur la mécanique de la rupture, il est
nécessaire de formuler un critère sur le choix de la direction. En fonction des conditions
de chargement (vitesse, température etc.), les mécanismes qui entrent en jeu au niveau de
la microstructure et le type de rupture sont différents (voir 1.1).

[Link] Comportement linéaire

Maximum de la contrainte d’ouverture En se basant sur les champs asymptotiques

présentés dans la partie [Link], Erdogan & Sih [ERD 63] proposent un critère simple
et intuitif qui postule que la fissure tend à se propager dans la direction qui maximise
le mode I. Ce critère est particulièrement adapté aux rupture fragiles par clivage. Ainsi,
si on se place dans le repère local cylindrique (F IG .1.17) et qu’on note σθθ la contrainte
circonférentielle, la propagation se fait dans la direction θc pour laquelle σθθ est maximale
à rayon fixé.

F IGURE 1.17: Repère local cylindrique en pointe de fissure.

Connaissant l’expression du champ de contrainte en pointe de fissure, il s’agit de

résoudre l’équation :
∂σθθ ∂2 σθθ
=0 et < 0. (1.20)
∂θ ∂θ2
Dans la mesure où la contrainte circonférentielle maximale est une contrainte principale,
la contrainte de cisaillement σrθ est nulle dans cette direction. On peut donc également
trouver θc en cherchant la solution de l’équation σrθ (θc ) = 0.

En reprenant l’expression de la contrainte de cisaillement plan asymptotique, il vient en

2D :
KI sin(θ) + KII (3cos(θ) − 1) = 0, (1.21)

Cette thèse est accessible à l'adresse : [Link]

ce qui après résolution donne :

  s 
2
1 KI KI 
θc = 2 arctan   − sign (KII ) 8 + . (1.22)
4 KII KII

Maximum du taux de restitution de l’énergie Les travaux de Erdogan & Sih

[ERD 63] étudient également la possibilité d’une propagation suivant la direction qui
maxime le taux de restitution de l’énergie G. Hussain [HUS 74] met en équation le calcul
de G pour une petite avancée de fissure δ orientée d’un angle θ par rapport à la fissure
existante (voir F IG .1.18). La relation obtenue est la suivante :
2 !θ
1 − πθ π

4 1
G(θ) = ∗
E 3 + cos2 θ 1 + πθ
× (1 + 3cos2 θ)KI2 + 8 sinθ cosθ KI KII + (9 − cos2 θ)KII2 , (1.23)

(
E
2 en déformations planes,
avec E ∗ = 1−ν
E en contraintes planes.

F IGURE 1.18: Calcul de G pour une avancée de δ et d’un angle θ.

L’angle qui maximise le taux de restitution de l’énergie est alors donné par :

∂G ∂2 G
=0 et < 0. (1.24)
∂θ ∂θ2

Minimum de la densité d’énergie de déformation On peut également citer les tra-

vaux de Sih [SIH 74] qui propose une direction de propagation suivant le minimum de la
densité d’énergie de déformation S.

[Link] Comportement non linéaire

Perte de stabilité matérielle Une autre approche proposée dans [BEC 02], [SÁN 08]
ou [HUE 09] consiste à supposer que la fissure se propage le long des plans de localisa-
tion de la déformation. En effet, à l’échelle locale, la perte de stabilité de l’état matériel

Cette thèse est accessible à l'adresse : [Link]

induit l’amorçage de la fissure. Cette direction est alors déduite par le biais d’une ana-
lyse de bifurcation. Pour identifier la direction de propagation il faut alors avoir accès à
l’opérateur tangent ce qui n’est pas toujours immédiat, et étudier la perte d’ellipticité de
cet opérateur. Ce critère se base en général sur l’utilisation d’une loi adoucissante comme
dans [LIN 09] (loi élasto-plastique avec adoucissement exponentiel) ou dans [CRE 13]
(loi élasto-plastique endommageable). Il est plutôt adapté à la propagation de fissure dans
des matériaux ductiles.

Maximum de la contrainte de cisaillement Récemment, Haboussa et al. ont déve-

loppé un critère bidimensionnel [HAB 12b] qui stipule qu’en présence d’une plasticité
non négligeable en pointe de fissure, la direction de propagation est pilotée par le maxi-
mum de la contrainte de cisaillement. En reprenant l’idée originale de facteurs d’inten-
sité de contraintes normalisées introduite par Schölmann et al. [SCH 02], ils ont ensuite
étendu le critère au cas tridimensionnel [HAB 12a] et choisissent cette fois-ci la direction
du maximum de la contrainte de von Mises.

1.2.4 Vitesse de propagation

D’après le critère développé par Erdogan and Sih (voir partie [Link]), le phénomène
de rupture est supposé piloté par l’intensité de la contrainte circonférentielle au voisinage
du front de fissure, évaluée en définissant le facteur d’intensité des contraintes circonfé-
rentielles Kθθ :
√
Kθθ = lim 2πr σθθ . (1.25)
r→0
Ce FIC est maximal pour l’angle θc présenté équation 1.22.

Tant que Kθθ reste inférieur à une valeur seuil KId appelée ténacité dynamique d’ini-
tiation, la rupture ne peut pas s’initier. Lorsque ce critère est atteint, la direction du maxi-
mum des contraintes circonférentielles θc définit la direction effective de propagation. Le
critère d’initiation s’écrit :
Kθθ < KId pas d’initiation, (1.26)
Kθθ = KId , θ = θc initiation. (1.27)
Après l’initiation, une autre valeur critique notée K1D (ȧ) est utilisée. Il s’agit de la
ténacité dynamique de propagation au dessus de laquelle une avancée de fissure est pro-
duite. Cette propagation reste effective tant que Kθθ est supérieur à la ténacité dynamique
d’arrêt notée KIa . Le critère de propagation peut être exprimé de la manière suivante :
Kθθ > KIa , Kθθ = K1D (ȧ) propagation. (1.28)
Ces différentes valeurs critiques sont abordées et caractérisées expérimentalement
par Grégoire et al. [GRE 07].

Cette thèse est accessible à l'adresse : [Link]

Il faut alors choisir une loi d’évolution de la propagation, c’est-à-dire une expression
reliant ce facteur d’intensité des contraintes (ou le taux de restitution de l’énergie) à la
vitesse d’avancée de la fissure ȧ. Les nombreuses expériences réalisées par le passé ont
montré qu’il existe une vitesse limite d’avancée du front de fissure et que cette vitesse
croît avec Kθθ . Ainsi Kanninen et Popelar [KAN 85] proposent une expression du type :
KIc
Kθθ = , (1.29)
1 − cȧR

où cR est la célérité des ondes de Rayleigh et KIc un paramètre matériau à déterminer

empiriquement.

La vitesse de propagation de la fissure s’écrit donc simplement :

KIc
ȧ = 1 − cR . (1.30)
Kθθ

Synthèse
Différentes approches de modélisation de la rupture ont été présentées dans cette par-
tie, certaines permettant d’appréhender une rupture en traction ou en cisaillement. On
pourra observer que ces travaux, bien que proposant de modéliser les deux types de rup-
ture, ne permettent pas le passage automatique de l’un à l’autre en cas de propagation.
C’est pour cette raison que nous avons fait le choix dans cette thèse de retenir l’approche
de Haboussa et al. qui bien que simplifiée permet une telle bascule.

Cette thèse est accessible à l'adresse : [Link]

1.3 Modélisation éléments finis d’une structure mince

Cette deuxième partie de l’étude bibliographique porte sur les aspects généraux
concernant les éléments finis de coques. Après avoir rappelé succinctement les étapes
de la résolution par éléments finis et introduit quelques notations, les hypothèses, la ci-
nématique et les précautions relatives à l’utilisation des éléments finis de coques sont
présentés.

1.3.1 Discrétisation par la méthode des éléments finis

Cette partie n’a pas pour objectif de rappeler la théorie des éléments finis mais de
fournir des éléments pour une meilleure compréhension des paragraphes suivants. Elle
permet également d’avoir les idées claires sur les étapes de construction d’un élément
fini et par la même occasion introduit quelques notations. Une présentation détaillée de
la méthode des éléments finis est faite dans [DHA 05], [BAT 90b], [BEL 99] ou encore
[ZIE 77].

La méthode des éléments finis est utilisée pour résoudre numériquement des équa-
tions aux dérivées partielles. Il s’agit donc d’une résolution approchée d’un problème,
où, grâce à la formulation variationnelle (forme intégrale des équations aux dérivées par-
tielles), les solutions du problème vérifient des conditions d’existence plus faibles que
celles des solutions du problème de départ et où une discrétisation permet de trouver une
solution approchée. Pour un problèmes de mécanique des structures, cette formulation
faible s’exprime à travers ce qu’on appelle le Principe des Travaux Virtuels (PTV) qui fait
intervenir quatre grandeurs : contraintes, déformations, forces et déplacements. Sa forme
intégrale s’exprime de la façon suivante :
Z Z Z Z
∀δ~u ¨ udV = −
ρ~u.δ~ σ̄¯ : δε̄¯ dV + f .δ~udV + T.δ~udS, (1.31)
D D D δD

avec :
δ~u un champ cinématiquement admissible à 0 ;
δε̄¯ le champ de déformation associé à δ~u ;
ρ la masse volumique ;
σ̄¯ le tenseur des contraintes de Cauchy ;
ε̄¯ le tenseur des déformations ;
f les forces volumiques ;
T les forces surfaciques.

Grâce à l’approximation nodale des déplacements, on passe de la recherche d’une va-

riable continue (le déplacement en tout point) à la recherche de variables discrètes (les

Cette thèse est accessible à l'adresse : [Link]

inconnues aux nœuds). Le lien entre les deux se fait via l’utilisation de fonctions d’inter-
polation le plus souvent polynomiales.

{~u(M)} = [N(M)] {Ue }. (1.32)

Connaissant le champ cinématique en tout point, il est donc possible d’accéder au champ
de déformation par dérivation, puis au champ de contrainte par la relation de comporte-
ment.

{ε(M)} = [B(M)] {Ue }, (1.33)

{σ(M)} = [H(M)] {ε(M)}, (1.34)

avec [B] la matrice d’opérateurs différentiels appliqués aux fonctions d’interpolation et

[H] la matrice de comportement liant contraintes et déformations.

Chacun des termes du PTV peut alors s’exprimer sous forme matricielle sur chaque
élément. On a alors : Z
¨ udV = {δUe }T [Me ]{Üe },
ρ~u.δ~ (1.35)
De

Z
avec [Me ] = [N(M)]T ρ[N(M)]dV matrice de masse élémentaire.
De

Z
σ̄¯ : δε̄¯ dV = {δUe }T [Ke ]{Üe }, (1.36)
De

Z
avec [Ke ] = [B(M)]T [H(M)][B(M)]dV matrice de raideur élémentaire.
De
Z Z
f .δ~udV + T.δ~udS = {δUe }T {Fe }, (1.37)
D δD

Z Z
T
avec {Fe } = [N(M)] { f }dV + [N(M)]T {T }dV .
De δDe

En remplaçant les équations précédentes dans l’équation 2.6, il vient :

[Me ]{ün } + [Ke ]{un } = {Fe } . (1.38)

Quelques techniques numérique sont nécessaires pour rendre systématique le calcul

de ces matrices. Dans un premier temps, plutôt que d’effectuer l’intégration directement
sur les éléments réels de la structure qui peuvent tous être différents en forme, rendant
l’intégration difficile, celle-ci est effectuée sur ce qu’on appelle un élément de référence.

Cette thèse est accessible à l'adresse : [Link]

Ce changement de domaine d’intégration se fait au moyen de la transformation géomé-

trique qui définit les coordonnées de tout point de l’élément réel à partir des coordonnées
du point correspondant de l’élément de référence.
Z Z
f (x, y, z)dxdydz = f (ξ, η, ζ)det[J]dξdηdζ, (1.39)
De Dref

avec [J] matrice de transformation géométrique, encore appelée matrice jacobienne.

Dans un second temps, il est nécessaire d’évaluer numériquement l’intégrale sur

l’élément de référence en la transformant en une somme discrète pondérée.
Z N pi
f dV = ∑ f (ξi )ωi . (1.40)
Dref i=1

F IGURE 1.19: Procédure de calcul d’un élément.

Pour résumer, la méthode des éléments finis revient à résoudre le PTV sur un ensemble
de sous-domaines et s’articule autour des étapes suivantes :
1. Discrétisation du milieu continu en sous domaines.
2. Construction de l’approximation nodale par sous domaine.

Cette thèse est accessible à l'adresse : [Link]

3. Calcul des matrices élémentaires.

4. Assemblage des matrices élémentaires.
5. Prise en compte des conditions aux limites.
6. Résolution du système d’équations.
Une vue synthétique des étapes 2 à 4 est présentée figure 1.19.

Maintenant qu’un bref aperçu des éléments finis a été réalisé, il convient de présenter
plus spécifiquement les particularités que présentent les éléments de type coque.

1.3.2 Cadre d’utilisation des coques

[Link] Généralités
La plupart des logiciels de simulation proposent des éléments coques, très utiles pour
modéliser des structures tridimensionnelles dont l’épaisseur est mince par rapport aux
autres grandeurs caractéristiques. On les retrouve dans différentes applications telles que
les carrosseries de voiture, les balles de tennis de table, les coquilles d’oeuf ou encore les
coques de bateau (voir F IG .1.20).

F IGURE 1.20: Modélisation de la coque d’un bateau.

L’idée de base dans la théorie de ces modèles de plaques et coques est d’utiliser
des hypothèses et simplifications physiquement justifiables à travers l’épaisseur pour
obtenir la déformation d’une structure tridimensionnelle élancée à partir d’un problème
formulé sur sa surface moyenne. On ramène alors ces problèmes tridimensionnels à une
étude bidimensionnelle. Un élément coque peut généralement ne pas être plan, saufs

Cette thèse est accessible à l'adresse : [Link]

pour certains éléments particuliers, dit éléments plaques. Dans la littérature, on trouve
profusion d’éléments finis de coques en passant des coques minces au coques épaisses.
Rien que dans le code industriel LS-Dyna, une quinzaine de formulations sont possibles.
Suivant l’orde de grandeur de h par rapport aux autres dimensions, on introduit souvent
l’adjectif "mince" ou "épais" aux coques. Ce qualificatif n’implique pas seulement
une caractéristique géométrique mais sous-entend également un rôle particulier des
déformations dites de cisaillement transverse (notée CT).

Les avantages à utiliser ce genre d’élément sont multiples. Ils permettent en premier
lieu de réduire de manière significative la taille des modèles numériques en terme de
nombre de degrés de liberté. De plus, quand on sait que lors d’un calcul en explicite
le pas de temps est piloté par la plus petite taille de maille (voir partie [Link]), une
discrétisation de l’épaisseur à l’aide d’éléments volumiques pénalise grandement le
temps de calcul. Ainsi, avec des éléments coques, l’épaisseur n’est plus directement
représentée et ne représente donc plus le facteur limitant pour le pas de temps critique, ce
qui conduit à réduire le temps de calcul.

En conception, les géométries de structures navales sont volumiques afin d’éviter des
erreurs de gabariage 1 . Passer au calcul éléments finis avec des éléments coque demande
donc un travail supplémentaire. Il est alors nécessaire de passer au feuillet de référence
de la structure qui est peut être composé du feuillet moyen, de la peau supérieure ou
inférieure de la géométrie. Néanmoins l’assemblage des tôles se voit facilité puisqu’il n’y
a plus qu’un seul bord à relier, contre quatre en volumique. Bien souvent dans l’industrie,
la géométrie d’une coque de forme quelconque est modélisée par un ensemble de facettes
planes.

[Link] Géométrie

Une plaque est un solide défini par une surface de référence plane et par une épaisseur
(notée h(x, y)) petite par rapport aux autres dimensions (longueur et largeur) de la
structure à modéliser (voir F IG .1.21).

F IGURE 1.21: Modélisation d’une plaque.

1. Contruction du gabarit du navire

Cette thèse est accessible à l'adresse : [Link]

On attache à la surface de référence (souvent confondue à la surface moyenne) un

repère orthonormé local (e1 , e2 , e3 ) associé au plan tangent de la structure différent du
repère global (X,Y, Z). La position des points de la plaque est donnée par les coordonnées
cartésiennes (x, y) de la surface moyenne et l’élévation z par rapport à cette surface,
comprise entre − h2 et h2 . Pour les coques, la position des points dépend des coordonnées
curvilignes notées x̃,ỹ et z̃.

On fait l’hypothèse des sections droites (non gauchissement), c’est à dire qu’un brin
normal au plan moyen reste droit et de longueur constante après déformation, ce qui
entraine notamment que la déformation dans l’épaisseur εzz est nulle. En d’autres termes,
on autorise le cisaillement transverse tout en négligeant la striction.

[Link] Cinématique
Pour la théorie des coques planes en petits déplacements et petites déformations, le
champ de déplacement UM d’un point M dans l’épaisseur de la plaque est donné dans le
repère local (e1 , e2 , e3 ) à partir des déplacements UG d’un point G de sa surface moyenne,
et la rotation θ = θx e1 + θy e2 la rotation de cette surface par rapport aux deux axes locaux
e1 et e2 (F IG .1.22). Ainsi on a :

UM = UG + θ ∧ z e3 = UG (x, y) + z θx (x, y)e2 (x, y) − z θy (x, y)e1 (x, y), (1.41)

qu’on peut également écrire :

     
UMx (x, y, z) UGx (x, y) −θy (x, y)
 UMy (x, y, z)  =  UGy (x, y)  + z  θx (x, y)  . (1.42)
UMz (x, y, z) UGz (x, y) 0

Pour faciliter la notation, on note aussi parfois la cinématique de la façon suivante :

     
UMx (x, y, z) UGx (x, y) βx (x, y)
 UMy (x, y, z)  =  UGy (x, y)  + z  βy (x, y)  , (1.43)
UMz (x, y, z) UGz (x, y) 0

avec βx = −θy et βy = θx .

Dans l’hypothèse des petites perturbations, le tenseur des déformations est défini par :

1
εi j = (ui, j + u j,i ). (1.44)
2
Ainsi, avec la cinématique introduite précédemment, les déformations tridimension-
nelles sont données en tout point par :

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.22: Cinématique d’un brin initialement perpendiculaire au plan moyen.



 εxx = ux,x + zβx,x ,
= uy,y + zβy,y ,

 εyy

γxy = (ux,y + uy,x ) + z(βx,y + βy,x ), (1.45)
= uz,x + βx ,



 γxz
= uz,y + βy ,

γyz

où γxy = 2εxy , γxz = 2εxz et γyz = 2εyz .

En décomposant la partie membrane de la partie flexion, on retrouve parfois la nota-

tion matricielle suivante :
{ε} = {e} + z{χ} + {γ}, (1.46)
avec
hεi = hεxx εyy γxy γxz γyz i,
hei = hux,x uy,y ux,y + uy,x 0 0i,
(1.47)
hχi = hβx,x βy,y βx,y + βy,x 0 0i,
hγi = h0 0 0 uz,x + βx uz,y + βy i,

hei : déformations de membrane ;

zhχi : déformations de flexion ; hχi : courbures ;
hγi : déformations de cisaillement transversal.

On observe qu’avec cette approche, les déformations membranaires (εxx , εyy et εxy )
varient linéairement dans l’épaisseur de la plaque, alors que les déformations de cisaille-
ment transverse (εxz et εyz ) sont elles constantes dans l’épaisseur.

Cette thèse est accessible à l'adresse : [Link]

Dans le cas courant où le matériau est élastique isotrope, nous avons :

 
1 ν 0
E  1 0 E
[H] = 1 0  ; [G] = G ; G= , (1.51)
1 − ν2 sym 1−ν 0 1 2(1 + ν)
2

où E et ν sont le module d’élasticité (module d’Young) et le coefficient de Poisson.

Pour résoudre le problème d’équilibre de la structure par la méthode des éléments

finis, les champs sont cherchés sous la forme donnée dans l’équation 1.32 :

UG = ∑ NA (ξ, η)UGA , θx = ∑ NA (ξ, η)θxA , θy = ∑ NA (ξ, η)θyA , (1.52)

A∈N A∈N A∈N

où G est un point de la surface moyenne et N l’ensemble des nœuds de l’élément.

[Link] Implémentation des coques courbes

Dans le cas d’une coque quelconque (non plane), on retrouve les même relations à la
différence près que ce sont les coordonnées curvilignes x̃, ỹ, z̃ qui interviennent :

UM = UG + θ ∧ z̃ e3 = UG (x̃, ỹ) + z̃ θx (x̃, ỹ)e2 (x̃, ỹ) − z̃ θy (x̃, ỹ)e1 (x̃, ỹ), (1.53)

Cette thèse est accessible à l'adresse : [Link]

ou encore :
     
UMx (x̃, ỹ, z̃) UGx (x̃, ỹ) −θy (x̃, ỹ)
 UMy (x̃, ỹ, z̃)  =  UGy (x̃, ỹ)  + z̃  θx (x̃, ỹ)  . (1.54)
UMz (x̃, ỹ, z̃) UGz (x̃, ỹ) 0
Pour l’implémentation des coques, il existe deux manières de procéder :

– La description curviligne des éléments coques : l’idée est de se placer dans le re-
père local de la surface moyenne (e1 , e2 , e3 ), puis on cherche les composantes cur-
vilignes des déplacements et rotations en tout point de coordonnée curviligne (x̃,ỹ).
Il est alors nécessaire de dériver les vecteurs dans la base locale si bien que les
courbures apparaissent explicitement. Cette description est très précise mais relati-
vement complexe à mettre en œuvre. Le champ de déplacement solution de l’équa-
tion 1.42 est recherché sous la forme :

UM = ∑ (NA (ξ, η)UGA ) + z̃ ∑ (NA (ξ, η)θxA ) e2 − z̃ ∑ NA (ξ, η)θyA e1 ,
A∈N A∈N A∈N
(1.55)
où e1 et e2 correspondent aux vecteurs tangents à la surface moyenne au point G.

– La description coques tridimensionnelles dégénérées ou massifs dégénérés : la

stratégie est de discrétiser non seulement la surface moyenne mais également la
normale à cette surface moyenne. Cela implique en particulier que le repère local
est obtenu en tout point par interpolation des repères locaux au nœuds de l’élément.
Le champ de déplacement solution de l’équation 1.42 est alors recherché sous la
forme :

UM = ∑ NA (ξ, η)UGA + z̃NA (ξ, η)θxA e2A − z̃NA (ξ, η)θyA e1A , (1.56)
A∈N
où e1A et e2A correspondent aux vecteurs tangents à la surface moyenne aux
nœuds.

1.3.3 Théorie de Kirchhoff-Love

Pour les coques dites minces, c’est à dire pour un rapport Lh > 20 avec L une longueur
caractéristique de la structure et h son épaisseur, les déformations dites de "cisaillement
transverse" γxz et γyz sont supposées nulles. On retrouve les mêmes hypothèses dans
la théorie des poutres d’Euler-Bernoulli adaptée aux poutres élancées. En reprenant
l’équation (1.45), on obtient :


 εxx = ux,x + zβx,x ,
 εyy = uy,y + zβy,y ,


γxy = (ux,y + uy,x ) + z(βx,y + βy,x ), (1.57)
0 = uz,x + βx ,




0 = uz,y + βy .


Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.23: Rotation d’une section droite initialement perpendiculaire à la

surface moyenne, théorie de Kirchhoff-Love

Cela suppose donc qu’une section orthogonale au plan moyen reste orthogonale après
déformation dans la configuration déformée : il s’agit de la théorie de Kirchhoff-Love
(F IG .1.23).

γxz = uz,x + βx = 0 βx (x, y) = −uz,x
⇒ (1.58)
γyz = uz,y + βy = 0 βy (x, y) = −uz,y

Ainsi les deux rotations de la normale à la surface moyenne sont liées aux déplace-
ments de la surface moyenne. Afin de garder une méthode conforme au sens mathéma-
tique, le modèle de Kirchhoff-Love doit être discrétisé par des éléments finis C1 . Ce type
d’élément fini est considéré a priori comme coûteux.

1.3.4 Théorie de Mindlin-Reissner

Pour les coques dites épaisses, c’est à dire pour un rapport Lh 6 10, le cisaillement
n’est plus négligeable et doit être pris en compte dans la cinématique. Cela se traduit
géométriquement par le fait qu’un brin normal au plan moyen peut ne pas lui rester ortho-
gonal après déformation comme on peut le voir figure 1.24. Ainsi, les rotations de la nor-
male, notées θx et θy sont des inconnues supplémentaires du problème (dans la théorie de
Kirchhoff-Love, ces rotations dépendent directement des autres inconnues cinématiques).
Dans ce cas, on peut faire l’analogie avec la théorie de Timoshenko pour les poutres avec
cisaillement transverse. Pour rappel, les déformations sont données par :

Cette thèse est accessible à l'adresse : [Link]



 εxx = ux,x + zβx,x ,
= uy,y + zβy,y ,

 εyy

γxy = (ux,y + uy,x ) + z(βx,y + βy,x ), (1.59)
= uz,x + βx ,

γ
 xz



γyz = uz,y + βy ,

avec βx = −θy et βy = θx .

F IGURE 1.24: Rotation d’une section droite initialement perpendiculaire à la

surface moyenne, théorie de Mindlin-Reissner.

En se référant aux expressions (1.48) et (1.59), il apparaît que les contraintes de ci-
saillement transverse (σxz et σyz ) sont supposées constantes dans l’épaisseur. Cela résulte
de la non dépendance des déformations de cisaillement transverse avec l’épaisseur. Ce-
pendant, du fait de cette formulation, il y a violation des conditions aux limites sur les
bords supérieurs et inférieurs. On a ici σxz 6= 0 et σyz 6= 0 au lieu de σxz = σyz = 0 . En
réalité, le champ de contraintes de cisaillement transverse est quadratique (voir F IG .1.25).
Afin de retrouver l’énergie interne du modèle que l’on aurait après résolution des équa-
tions d’équilibre tridimensionnel avec le champ analytique, on fait apparaître un facteur
de correction de cisaillement k entre l’énergie interne de cisaillement transverse analy-
tique et celle donnée par la théorie de Mindlin :

Wxz(Analytique) = k ×Wxz(Mindlin) . (1.60)

Ce facteur de correction est pris égal à 56 . On peut montrer qu’ainsi l’énergie de

cisaillement transverse est la même avec les deux représentation du cisaillement dans
l’épaisseur. Récemment, Bouclier [BOU 12] a montré de manière générale comment

Cette thèse est accessible à l'adresse : [Link]

(a) (b)

F IGURE 1.25: Forme du champ de cisaillement transverse dans l’épaisseur dans le

cas analytique (a) et dans la théorie de Mindlin (b).

on pouvait formuler directement une coque sans verrouillage à partir d’une formulation
mixte.

Le modèle de Mindlin-Reissner ne nécessite que des éléments finis C0 pour assurer

la conformité mathématique. Ils sont donc moins "couteux" que les éléments du type
Kirchhoff-Love et la plupart des éléments coques proposés dans les logiciels du com-
merce reposent sur cette théorie. Néanmoins, ces éléments sont sujet au phénomène dit
de "verrouillage numérique" qui vient fausser les résultats.

1.3.5 Verrouillage numérique

L’expérience montre que la discrétisation par éléments finis surestime souvent la
rigidité d’un modèle : on obtient des déplacements plus faibles que ceux qu’aurait prédit
la mécanique des milieux continus. Cet écart est généralement léger, mais il arrive
toutefois que certains éléments présentent une rigidité extrêmement élevée, supérieure
de plusieurs ordres de grandeurs à ce qu’elle devrait être. Ce phénomène s’appelle le
verrouillage et apparaît souvent pour les plaques dont le rapport de l’épaisseur sur la
longueur est faible. Des exemples de verrouillage sont présentés dans [BAT 90a] ou
encore [LAS 09].

Remarque : La simulation du comportement des matériaux incompressibles, comme les

polymères, peut elle aussi faire apparaître des phénomènes de verrouillages.

Origine du verrouillage En coque, le verrouillage se manifeste d’autant plus que la

structure est sujette à de la flexion. En effet, les fonctions d’interpolation peuvent ne pas
être capables de reproduire correctement la flexion sans cisaillement, ou la flexion sans
membrane des coques minces. Dans ce cas, des contraintes supplémentaires résultat de
termes erronés en cisaillement et membrane apparaissent et absorbent la majeure partie
de l’énergie de déformation alors que ce sont justement ces contributions qui devraient

Cette thèse est accessible à l'adresse : [Link]

être quasiment nulles. Pour résumer, la solution aura plutôt tendance à se bloquer sur
la solution de cisaillement ou de membrane plutôt que de décrire le comportement en
flexion.

Le plus violent est le verrouillage en cisaillement transverse qui apparait déjà en

plaque et est indépendant des conditions aux limites. Ensuite, vient le verrouillage en
membrane qui est fonction des conditions aux limites puisqu’il apparait lorsque l’on mo-
délise un problème de type flexion pure. La matrice de rigidité d’un élément fini de coque
peut être considérée comme la somme des matrices associées aux termes de membrane,
de flexion, de couplage membrane/flexion (nulle sous certaines hypothèses) et de cisaille-
ment transverse. On convient de noter son expression, en faisant apparaitre la dépendance
d’ordre un en h, sous la forme :
[K] = h[Km ] + h[Kc ] + h2 [Km f ] + h3 [K f ]. (1.61)

Verrouillage en cisaillement transverse Pour le cas du cisaillement transverse, on

considère un problème de plaque en flexion et le modèle de Mindlin/Reissner. Ainsi,
l’équation 1.61 devient :
[K] = h[Kc ] + h3 [K f ]. (1.62)
Il est alors nécessaire que [Kc] tende vers 0 quand h tend vers 0 afin d’éviter que
le terme cisaillement transverse devienne prépondérant dans l’expression de la rigidité.
Dans le cas contraire, le verrouillage apparait du fait de la dépendance en h des matrices.
L’élément conduit ainsi à une solution influencée uniquement par le cisaillement au
lieu d’être associée au phénomène de flexion Il y a dans ce cas là moins d’énergie pour
déformer la plaque, ce qui donne une solution erronée trop rigide. C’est l’apparition de
cette rigidité artificielle qui est appelée "verrouillage en cisaillement" (shear locking).
Cela apparait de manière significative pour les plaques qui ont un rapport Lh > 20.

L’utilisation d’une technique isoparamétrique ne permet pas aux déformations de ci-

saillement transverse de devenir négligeables (voir équation 1.59) puisque la flèche et
les rotations sont approximées dans le même espace polynomial ; la somme d’un terme
constant et d’un terme linéaire ne peut être nul en tout points.

Verrouillage en membrane En coque et en considérant un problème de flexion pure,

la même constatation peut être faite. On considère alors le modèle Kirchhoff Love négli-
geant les déformations de cisaillement transverse. L’équation 1.61 devient alors :
[K] = h[Km ] + h2 [Km f ] + h3 [K f ]. (1.63)
Si les déformations de membrane ne sont pas négligeables, alors le comportement asymp-
totique de [K] sera dominé par les termes de membrane alors que c’est un problème de
flexion. On aura alors verrouillage de part la dépendance en h.
Des éléments qui minimisent ce verrouillage ont été développés, et permettent de
traiter ces phénomènes.

Cette thèse est accessible à l'adresse : [Link]

Traitement du verrouillage en cisaillement transverse Afin de remédier au ver-

rouillage numérique des plaques élancées dû au cisaillement, différentes techniques ont
été proposées :
1. L’intégration réduite sélective (SRI) : La méthode consiste à combiner astucieuse-
ment les techniques d’intégration complète et réduite (matrice de rigidité de flexion
et membrane en intégration complète et matrice de rigidité en cisaillement en inté-
gration réduite).
2. L’intégration réduite : La méthode consiste à prendre un seul point d’intégra-
tion de Gauss (pour cet élément) au centre du plan moyen au lieu de quatre points
comme en intégration complète. L’avantage se résume principalement à un gain en
temps de calcul tout en levant le verrouillage. Néanmoins, les éléments à intégra-
tion réduite peuvent conduire à développer des modes à énergie nulle bien qu’il y
ait des déplacement aux noeuds. Cela peut entrainer de fortes distorsions, créant des
zones où aucune contrainte n’accompagne la déformation constatée de l’élément.
On parle alors du phénomène d’Hourglass (ou sablier) caractérisé par les formes
présentées dans la figure 1.26. Plusieurs méthodes de stabilisations ont été intro-
duites par Belytschko et al. [BEL 99] afin d’empêcher l’apparition de ces modes
d’Hourglass.

F IGURE 1.26: Modes de Hourglass associés à une énergie de déformation nulle.

3. Champ assumé de déformation de cisaillement transverse : C’est une méthode

qui peut remplacer les deux précédentes. La particularité de celle-ci est de calculer
le cisaillement transversal γ d’un quadrangle non pas directement aux quatre points
d’intégration comme pour la membrane et la flexion mais de les interpoler à l’aide
des valeurs calculées aux quatre points milieux des côtés du quadrilatère. Les
composantes suivant xz et yz servent ensuite à définir un champ de déformations
"assumé" dans les directions ξ et η en un point courant. Cela permet d’éviter le
phénomène de verrouillage en cisaillement transverse.

Cette liste n’est pas exhaustive et on peut par exemple rajouter l’utilisation d’une
formulation mixte permettant de diminuer l’influence prépondérante du cisaillement par
rapport à la flexion.

Cette thèse est accessible à l'adresse : [Link]

1.3.6 Quelques éléments finis de coques classiques

Eléments de Kirchhoff Ces éléments finis sont basés sur la théorie des plaques minces
en flexion, appelée théorie "continue" de Kirchhoff qui ne prend pas en compte les ef-
fets de cisaillement transverse. Un élément conforme au sens mathématique exige une
continuité C1 , c’est à dire la continuité de Uz et de Uz,n le long du contour des éléments.
Il n’est cependant pas simple de satisfaire cette continuité avec un champ Uz défini sous
forme polynomiale sur un élément n’ayant comme variables nodales que Uz et ses déri-
vées premières (pour représenter la flexion). Dans la littérature, on peut citer l’élément
rectangulaire R16 de type Hermite, les éléments triangulaires HCT , T 18, HSM6 ou en-
core HSM9 qui sont présentés dans [BAT 90a].

Éléments de Kirchhoff Discret Plusieurs éléments finis permettent d’obtenir les so-
lutions numériques des problèmes de plaques minces en accord avec les hypothèses de
Kirchhoff, mais cette fois-ci en s’affranchissant de la continuité C1 . On définit donc des
éléments dit de type "Kirchhoff-discret" à faible nombre de noeuds et de degrés de li-
berté sans utiliser la théorie continue de Kirchhoff. Les hypothèses de Kirchhoff sont en
fait introduites sous forme particulière sur le contour, ou éventuellement à l’intérieur de
l’élément. Le triangle DKT (Discret Kirchhoff Triangle) et le quadrangle DKQ (Discret
Kirchhoff Quadrangle) sont construits sur ce principe (voir [BAT 90a]).

Éléments de Mindlin Reissner Ces éléments finis sont basés sur la théorie de Mindlin
Reissner et sont valables à la fois sur les plaques fines et épaisses. Comme énoncé pré-
cédemment, ces éléments sont sujets au phénomène de verrouillage numérique qu’il faut
traiter afin d’obtenir des résultats cohérents. L’élément MITC4 (Mixed Interpolation of
Tensorial Components) [DVO 84, BAT 86] l’élément S4R d’Abaqus [SYS 12] ou encore
l’élément Q4γ sont des éléments basés sur ces hypothèses.

Éléments de type massif coque L’idée de base de ces éléments est de proposer une
formulation massive pour modéliser une structure mince. Il s’agit d’éléments tridimen-
sionnels à intégration réduite sauf dans une direction privilégiée appelée épaisseur. La
formulation est plus simple que dans le cas d’un élément de coque puisqu’elle n’utilise
que des degrés de liberté de déplacement dans un repère cartésien. il est également très
facile de les coupler avec des éléments volumiques classiques et la loi de comporte-
ment n’est plus écrite en contraintes planes. On peut citer la famille d’éléments SHB
[ABE 02, LEG 02].

Cette thèse est accessible à l'adresse : [Link]

1.4 Modélisation numérique de la propagation de fissure

Cette partie est consacrée aux différents outils présent dans la littérature et dont dis-
pose le numéricien pour modéliser la propagation de fissure. Elle donne un aperçu des
principales méthodes numériques utilisées pour modéliser ces phénomènes.

1.4.1 Méthodes non basées sur les éléments finis

[Link] Méthode des éléments de frontières
Cette méthode exposée dans [ALI 02, YAN 95, TUH 97] présente comme principal
avantage de ne demander qu’une discrétisation de la frontière du domaine et non de l’in-
térieur (voir F IG .1.27). Elle permet une description aisée de la fissure et son évolution est
naturelle. Les équations aux dérivées partielles sont transformées en des équations inté-
grales sur la frontière par le biais d’une solution fondamentale dite de Green. Cette mé-
thode n’est applicable que si une telle solution existe ce qui n’est pas le cas des problèmes
élastoplastiques. Pour ces problèmes, un maillage de la zone plastifiée est nécessaire ce
qui revient à renoncer au principal avantage de cette méthode. Davantage d’informations
concernant les problèmes liés à cette méthode sont précisées dans [POR 92] et [FED 04].

F IGURE 1.27: Trajet de fissure pour différentes directions de charge imposée

[FED 04].

[Link] Méthodes sans maillage

Parmi les méthodes sans maillage, la méthode SPH (Smoothed Particule Hydro-
dynamics) [GIN 77] [LUC 77] et la méthode "Element Free Galerkin" [BEL 94] sont
certainement les plus utilisées. La première se base sur une formulation forte évaluée
en un certain nombre de points et la deuxième se base sur une formulation faible de

Cette thèse est accessible à l'adresse : [Link]

type intégrale. Ces méthodes sont très efficaces lorsque l’on traite de très grandes
déformations ou de fragmentation. Elle ne raisonnent plus sur un découpage de la
géométrie en éléments, mais sur une représentation du volume par une densité de noeuds
dont la connectivité pourra évoluer au cours des calculs.

Alors qu’avec la méthode des éléments finis classiques seuls les noeuds de l’élément
où l’on cherche à évaluer le champ interviennent, ces méthodes se basent sur l’approxi-
mation d’un champ faisant intervenir tous les noeuds présent dans un voisinage (voir
F IG .1.28). L’influence de chacun de ces noeuds est pondérée par une fonction poids.
L’avantage principal de ces méthodes sans maillage est l’absence de discrétisation de la
fissure. Des développements en mécanique linéaire et non-linéaire de la rupture proposent
des résultats intéressants [XU 98], [FLE 97] et [DUF 04].

F IGURE 1.28: Représentation des noeuds et de leur domaine d’influence

Cependant, ces méthodes sont couteuses en temps de calcul ; d’une part du fait de
la recherche systématique de l’influence des noeuds voisins et d’autre part puisqu’une
discrétisation plus fine qu’avec les éléments finis classiques est nécessaire. De plus,
la méthode SPH classique souffre de quelques problèmes de stabilité. Le premier
concerne une instabilité en tension qui apparait lorsque les nœuds SPH s’éloignent trop
les uns des autres. La densité de points dans le support du noyau devient de plus en
plus faible et l’approximation devient plus grossière. Ajouté à cela, un problème de
sous-intégration se manifeste du fait que toutes les grandeurs cinématiques sont stockées
sur une même particule. Des modes comparables aux modes d’Hourglass en éléments
finis émergent conduisant à la formation d’agglomérats. La formation de ces agglomérats
rend la répartition de billes moins homogène et l’éloignement progressif de ces billes va
entraîner une baisse de l’intensité de leur interaction ce qui aura pour effet à son tour
d’augmenter l’éloignement. Les billes vont finir par s’éloigner suffisamment pour sortir
de leur voisinages respectifs et donc perdre totalement leur interactions. La barre va
alors se fracturer artificiellement comme on peut le voir figure 1.29. Enfin, les dernières
difficultés liées à ces méthodes sont dans l’application des des conditions aux limites

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.29: Fracture artificielle lors d’un essai de traction avec des élément SPH
classiques (tiré de [MAU 08b]).

cinématiques.

Pour résoudre le problème d’instabilité en tension, Belytschko et al. [BEL 00] pro-
posent une stratégie basée sur une formulation lagrangienne totale plutôt que la formu-
lation classique. Ce sont alors les positions initiales des billes qui interviennent dans les
calculs d’interaction et l’éloignement des particules sur la configuration déformée n’en-
gendre donc plus de problème d’instabilité. Pour ce qui est des phénomènes liés à la
sous-intégration (qui sont d’autant plus gênant pour les formulations coque SPH, voir
[MAU 08a]), la méthode des stress points introduite par Dyka [DYK 97] semble fournir
des résultats satisfaisant.

F IGURE 1.30: Comparaison des simulations SPHC et EF sur un essai de

perforation [CAL 11, RUS 09].

Certains travaux récents traitent d’une formulation de coques couplée à la méthode

SPH (F IG .1.30) pour des applications de fissuration [MAU 08a, MAU 08b, CAL 11].

1.4.2 Méthodes basées sur les éléments finis classiques

[Link] Méthode de remaillage
Une première manière de représenter une fissure dans un maillage élément finis est
tout simplement de l’introduire explicitement en maillant directement ses lèvres et sa
pointe (F IG .1.31). Ainsi, l’état initial correspond à un maillage de la structure compa-
tible avec la préfissure voulue. Ensuite, pour faire avancer cette fissure, les positions de la

Cette thèse est accessible à l'adresse : [Link]

pointe et des lèvres sont mises à jour ce qui implique une étape de remaillage de la struc-
ture. L’analyse numérique des champs mécaniques d’une pièce fissurée est fortement liée
à la qualité du maillage de cette pièce, et plus particulièrement du voisinage de la pointe
de fissure. Cela implique donc l’utilisation d’un maillage suffisamment fin au voisinage
de la pointe de fissure. Lors de l’étape de remaillage qui est souvent complexes à mettre
en œuvre, il est nécessaire de projeter les différents champs du pas précédent sur le nou-
veau maillage ce qui conduit à une non conservation de l’énergie du système. On peut par
exemple citer les travaux récent de Chiaruttini et al. [CHI 12] qui se basent cette méthode
de remaillage automatique.
Cette méthode semble discutable et n’est jamais utilisée en dynamique rapide. Elle est
peu employée dans l’industrie car elle est généralement trop coûteuse en temps de calcul.

F IGURE 1.31: Propagation de fissure par une technique de remaillage automatique

[BOU 03].

[Link] Méthode d’élimination d’éléments

F IGURE 1.32: Propagation d’une fissure par érosion.

La méthode d’élimination d’éléments, souvent appelée méthode d’érosion des élé-

ments, est l’un des outils numériques les plus simples utilisés pour simuler la propaga-
tion des fissures en se basant sur des éléments finis classiques. Cette méthode consiste

Cette thèse est accessible à l'adresse : [Link]

à retirer définitivement un élément lorsqu’un critère de rupture est atteint. Ce critère est
fréquemment piloté par un endommagement ou une déformation plastique cumulée. Ri-
goureusement, cette méthode viole le principe de conservation de la masse et de l’énergie
bien qu’en général les éléments érodés sont de petite taille et leur masse est minime. Pour
limiter ce problème, certains codes comme LS-Dyna proposent de remplacer les éléments
érodés par des particules SPH qui peuvent évoluer avec ou indépendamment de la struc-
ture.
Un autre inconvénient réside dans la représentation de la fissure. Comme on peut
l’observer sur la figure 1.32, celle-ci ne peut d’une part que suivre le maillage et d’autre
part, elle n’est pas d’épaisseur nulle : l’écartement des lèvres de la fissure est égale à la
taille de l’élément érodé.

F IGURE 1.33: Représentation de la fissure par des éléments érodés [SON 08].

De plus, la suppression d’un élément entraîne une libération brutale de l’énergie dis-
sipée plastiquement dans celui-ci, ce qui conduit parfois à une propagation artificielle de
la fissure. Ainsi, il est préférable d’ajouter une loi constitutive dans laquelle la contrainte
tend vers zéro pour des déformations suffisamment importantes. La fissure est alors mo-
délisée par un groupe d’éléments "désactivés" (voir F IG .1.33) qui sont le plus souvent
remplacés par des éléments à contrainte nulles, c’est-à-dire sans résistance mécanique.
Des exemples de telles lois sont montrées figure 1.34. Cependant, l’implantation de tels

F IGURE 1.34: Courbes d’évolution contraintes-déformations pour un matériau

endommageable montrant un adoucissement élastique et un écrouissage plastique
[SON 08].

Cette thèse est accessible à l'adresse : [Link]

comportements matériaux reste problématique en dynamique de la rupture. En effet, les

simulations numériques effectuées avec des matériaux adoucissant ou des matériaux en-
dommageantes donnent des résultats qui ne sont pas objectifs vis à vis du maillage. Cela
se traduit généralement par une localisation (voir [BEN 93, FOR 04]) artificielle de la
déformation ou de l’endommagement dans une zone dont la taille tend vers zéro lorsque
le maillage est raffiné [BAZ 76]. En particulier, la dissipation d’énergie par fissuration
y devient nulle, ce qui est physiquement irréaliste. Ainsi, l’orientation et la finesse du
maillage influencent de manière significative les résultats comme on peut le voir dans la
figure 1.35. La comparaison de quelques résultats expérimentaux avec cette méthode par
Song et al. [SON 08] a mis en doute sa capacité pour être appliquée au cas dynamique.

F IGURE 1.35: Localisation pathologique de l’endommagement, effet de la taille et

de l’orientation du maillage [FOR 04].

Afin de régulariser ce problème et ainsi d’obtenir une solution objective, c’est-à-dire

indépendante du maillage, plusieurs familles de limiteurs de localisation ont été intro-
duites. Une bonne synthèse de ces limiteurs a été récemment faite dans [MOR 14]. On
peut par exemple citer les modèles à conservation de l’énergie de fissuration qui per-
mettent de faire en sorte que l’énergie libérée par la création de la surface de la fissure
soit égale au taux de restitution d’énergie critique du matériau (voir paragraphe [Link]). Il
existe également des modèles dits non-locaux qui permettent de régulariser le problème en
introduisant une longueur caractéristique. Une première approche dite intégrale consiste

Cette thèse est accessible à l'adresse : [Link]

à remplacer une variable locale a par sa "moyenne" ã sur le voisinage du point considéré :
R
ΩRω(x, y)a(y)dy
ã = , (1.64)
Ω ω(x, y)dy

où ω est une fonction de pondération qui détermine l’influence du point y sur le point
d’intérêt x. Cette fonction de pondération est souvent une gaussienne. Cette approche qui
a été initialement développée pour des modèles d’endommagement par Pijaudier-Cabot
et Bazant [PIJ 87] est aujourd’hui également utilisée pour des matériaux élasto-plastiques
et élasto-plastiques-endommageables [BAZ 02]. L’approche différentielle est une autre
forme de modèle non-local et consiste à introduire le second gradient de la défor-
mation ou bien d’une variable interne dans les équations. Là encore, des travaux ont
été menés aussi bien sur des modèles élasto-plastiques [AIF 84, LAS 88, DEB 92]
qu’endommageable [PEE 96, FRÉ 96]. On trouve également les modèles à effet de
vitesse introduits par Needleman [NEE 88] qui propose d’ajouter un terme de viscosité
dans la loi de comportement. Cela a également pour effet d’améliorer la stabilité des
simulations numériques. Enfin, dans le même esprit, on trouve les modèles à effet retard
initialement proposé par Ladevèze et al. [LAD 92] puis complétés pour devenir des
modèles à taux limité [DEU 97, SUF 04]. En se basant sur le fait qu’une fissure ne peut
pas apparaître et se propager instantanément dans une structure, il introduit une viscosité
dans la loi d’évolution de l’endommagement. Cela a pour effet de limiter la croissance
de l’endommagement et donc celui-ci a le temps de se développer sur une largeur de
plusieurs éléments avant d’atteindre une valeur critique.

En résumé, ce type de méthode est très utilisée pour sa simplicité de mise en œuvre
et vise davantage à étudier le comportement macroscopique de la structure dans un cadre
conservateur. Dès lors que l’on souhaite connaître plus finement la réponse, l’utilisation
d’une autre méthode notamment plus objective est préférable.

[Link] Méthode des éléments cohésifs

Comme décrit dans le paragraphe [Link], les modèles de zones cohésives fournissent
un cadre sûr et robuste pour prédire la cinétique de propagation lorsque la localisation
géométrique des fissures potentielles est connue. Cette approche est intermédiaire entre
les modèles basés sur la mécanique de la rupture utilisant des critères d’amorçage (taux
de restitution d’énergie ou critères en contrainte) et ceux utilisant la mécanique continue
de l’endommagement.

Néanmoins, cette méthode possède quelques inconvénients. Tout d’abord elle modifie
la rigidité globale de la structure en introduisant un élément d’épaisseur nul mais de
rigidité finie. Cela perturbe entre autre le trajet des ondes mécaniques en dynamique. De
plus, il est nécessaire de connaître a priori le trajet de la fissure en disposant des éléments
cohésifs le long de ce trajet. Par ailleurs, comme il a été montré dans [XU 94, ZHO 04],
ce modèle n’est pas objectif vis à vis du maillage. La direction de propagation montre

Cette thèse est accessible à l'adresse : [Link]

alors une forte dépendance au maillage en terme d’orientation et de taille des éléments.

Une évolution de cette méthode est l’utilisation des zones cohésives dans le cadre
d’une méthode de partition de l’unité telle que la méthode X-FEM (développée dans le
paragraphe [Link]) comme c’est le cas dans [MOË 02] ou plus récemment pour le cas
spécifique de la déchirure ductile [SIM 10]. Cela permet entre autre à la fissure d’évoluer
librement au sein du maillage sans nécessairement suivre le pavage. Un critère d’insertion
des zones cohésives ainsi que qu’un critère d’orientation doivent alors être définis. Pour
la propagation, les critères sont basés sur des grandeurs internes à la zone cohésive. Afin
de calculer le comportement dans cette zone, des points de Gauss sont répartis sur la
surface de la fissure, puis les contraintes cohésives sont ensuite interpolées pour obtenir
la contribution de la zone cohésive aux efforts nodaux.

1.4.3 Méthodes basées sur la partition de l’unité

[Link] Généralités

L’idée sous-jacente de ces méthodes est d’améliorer la qualité de la solution fournie

par la méthode des éléments finis non pas en raffinant la discrétisation, mais en injectant
des fonctions a priori proches de la solution attendue dans la base des fonctions utilisée
pour évaluer le problème. Il est ensuite nécessaire de formuler un critère de propagation
se basant sur cette discrétisation enrichie.

Ce concept introduit dans la fin des années 90 par Melenk et Babuska [MEL 96] s’ap-
plique donc particulièrement bien à tout problème présentant des singularités, y compris
la mécanique de la rupture. Pour cette application, on peut citer la méthode des éléments
finis étendus X-FEM de Belytschko, Moës et al. [MOË 99, DOL 00a, BLA 99], la mé-
thode des éléments finis généralisés proposée par Strouboulis et al. [STR 00, DUA 00],
ainsi que des méthodes utilisant des zones cohésives développées par de Borst et al.
[BOR 04, WEL 01, REM 03]. En pratique et pour des raisons de coût de calculs, on
limite l’introduction de fonctions d’enrichissement où cela est nécessaire. Ainsi on crée
une couche d’éléments de transition, appelés "blending elements" et visibles figure
1.36, pour lesquels une partie seulement des noeuds est enrichie et qui ne vérifie donc
plus vraiment la propriété de partition de l’unité. La présence de cette zone et la façon
dont elle est traitée peut avoir une influence sur l’ordre de convergence de la méthode
[CHE 03].

Soit un domaine Ω maillé par un ensemble d’éléments E et discrétisé par un ensemble

de nœuds N . Un ensemble de fonctions de forme NA s’appuie sur ces ensembles. On peut
alors représenter une approximation éléments finis d’un champ u de la façon suivante :

u(x) = ∑ NA (x)uA . (1.65)

A∈N

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.36: Mise en évidence des "blending elements" sur des domaines Ω
discrétisés en 1D et 2D [FRI 08].

Mathématiquement, si les fonctions de forme NA constituent une partition de l’unité dans

le domaine Ω, on a :
∑ NA (x) = 1, ∀x ∈ Ω. (1.66)
A∈N

On peut alors enrichir l’approximation de u avec une fonction d’enrichissement φ(x) de

la façon suivante :
u(x) = ∑ NA (x)uA + ∑ NB (x)φ(x)uB , (1.67)
A∈N B∈N e

où N e représente l’ensemble des nœuds auxquels on choisit de placer des degrés de

libertés enrichis uB . La fonction d’enrichissement est choisie de manière à capturer
plus efficacement la solution du problème considéré. Selon le problème traité, il peut
parfois être judicieux d’introduire plusieurs fonctions d’enrichissement si une unique
fonction ne permet pas de capturer totalement la discontinuité. Cette modification
de la base des fonctions utilisées pour résoudre le problème permet d’introduire des
discontinuités fortes telle qu’une fissure dans les champs de déplacement et déformations.

Les principaux intérêts de ces méthodes sont l’amélioration de l’ordre de convergence

et le découplage entre la discontinuité et le maillage. En effet, la discontinuité étant conte-
nue dans les fonctions d’approximation, sa géométrie apparait de manière implicite. Ainsi,
pour la faire évoluer au cours du temps, l’étape de remaillage n’est plus nécessaire.

Cette thèse est accessible à l'adresse : [Link]

[Link] La méthode des éléments finis étendus X-FEM

En mécanique de la rupture, c’est la fissure et son effet sur les différents champs
(déplacement, déformation etc.) que l’on souhaite représenter. Ainsi, en se basant sur
la méthode de la partition de l’unité présentée au paragraphe précédent, les fonctions
pour enrichir la base des solutions doivent être capable de représenter d’une part la
discontinuité du champ de déplacement au passage de la fissure loin de la pointe, et
d’autre part la singularité du champ de déplacement autour de la pointe.

La méthode des éléments finis étendus (ou "eXtended Finite Element Method") pro-
posée par Belytschko et Black [BLA 99] est basée sur ce principe. Pour représenter les
lèvres de la fissure, des degrés de liberté représentant le saut de déplacement sont ajou-
tés [MOË 99]. Ces degrés de liberté sont obtenus en multipliant une fonction de forme
élément fini classique par une fonction de type Heaviside H introduite de la manière
suivante :
1 si x est au dessus de Γ
H (x) = (1.68)
−1 si x est en dessous de Γ
où Γ représente la géométrie de la fissure et x le vecteur position d’un point de Ω. L’équa-
tion 1.67 devient alors :

ui (x) = ∑ NA (x)uiA + ∑ NB (x)H (x)uiB , (1.69)

A∈N B∈N cut

où N est l’ensemble des nœuds du maillage considéré, NA (x) est la fonction de forme
éléments finis classique associée au noeud A, uiA est la ième composante du degré de
liberté associé, N cut est l’ensemble des noeuds appartenant à un élément traversé par la
fissure et uiB est la ième composante du degré de liberté enrichi associé au noeud B.

Sur le même principe, le champ singulier au voisinage de la pointe peut être approché
à l’aide d’une base de fonctions d’enrichissement représentant exactement les champs
asymptotiques du problème d’une fissure dans un milieu élastique. On enrichi alors l’élé-
ment contenant la pointe de fissure à l’aide de fonctions singulières (voir F IG .1.37). De
nombreux travaux [BLA 99, MOË 99, DUA 00] ont montré que la base √ la plus adéquate
était composée de quatre fonctions singulières avec une évolution en r :
√ n φ φ φ φ o
{Fk }k=1..4 = r sin( ), cos( ), sin( )sin(φ), cos( )sin(φ) , (1.70)
2 2 2 2
où (r,θ) sont les coordonnées dans le repère local lié front de la fissure, d’un point appar-
tenant à un plan normal au front. Le champ de déplacement est alors décrit par l’approxi-
mation éléments finis étendus suivante :

ui (x) = ∑ NA (x)uiA + ∑ NB (x)H (x)uiB + ∑f ront ∑ NC (x)Fk (x)uiC,k , (1.71)

A∈N B∈N cut C∈N k
| {z } | {z }
Standard Enrichissement

Cette thèse est accessible à l'adresse : [Link]

où N f ront est l’ensemble des nœuds d’éléments contenant le front de fissure et uiC,k la
ième composante du degré de liberté enrichi associé au noeud C et à la kème fonction
d’enrichissement.

F IGURE 1.37: Stratégie d’enrichissement pour une fissure quelconque placée sur
un maillage [ELG 06a].

Historiquement, ce sont les fonctions singulières qui ont été introduite les premières.
Les fonctions saut sont ensuite apparues afin de s’affranchir de la dépendance au maillage.
Actuellement la méthode X-FEM ne s’applique que sur des structure déjà fissurées.
L’avantage principal de cette méthode est de permettre à la fissure d’évoluer sans néces-
saire suivre le maillage. A l’heure actuelle, on ne peut pas parler de modèle indépendant
du maillage car les critères de propagation proposés dans la littérature (voir Partie 2.3.1)
nécessitent des maillages suffisamment fin au voisinage de la pointe de fissure. Ceci est
indispensable à la fois pour capter "suffisamment de physique" pour faire un choix de
direction de propagation par exemple. De plus, si on fait le choix de ne garder que l’en-
richissement de type Heaviside, la pointe de fissure se trouve nécessairement sur le bord
d’un élément et l’avancée ne peut pas être inférieure à la taille de l’élément suivant.

1.4.4 Méthodes diffuses

[Link] Méthode à champ de phase
Les difficultés à simuler numériquement la propagation de chemins de fissuration
complexes à l’aide de la mécanique de la rupture ont conduit certains auteurs à utiliser
des méthodes à champ de phase, souvent appelées méthodes phase field. Ces méthodes
ont initialement été introduites pour étudier l’évolution de la microstructure [WAR 95]

Cette thèse est accessible à l'adresse : [Link]

puis ont été reprises pour traiter la propagation dynamique de fissure [KAR 01, MIE 10].

Le principe de l’approche est d’introduire un champ auxiliaire similaire à un champ

d’endommagement. Le champ prend la valeur 1 sur la fissure et décroît régulièrement vers
0 pour des points matériels s’éloignant de la fissure. L’approximation régularise ainsi une
discontinuité géométrique dans le solide par des zones diffuses d’apparition de la fissure.
On peut rapprocher cette méthode des modèles à gradient.
On peut citer les travaux récents de Borden et al. [BOR 12] qui introduisent l’utilisa-
tion des phase field dans des éléments isogéométriques pour étudier la rupture fragile en
dynamique.

[Link] Méthode Thick Level Set

Récemment proposée par [MOË 11], une nouvelle manière de faire évoluer l’en-
dommagement dans les solides a vu le jour. La méthode des level set épaisses, encore
appelée Thick Level Set ou TLS, a été formulée afin d’éviter la localisation non physique
observée avec les modèles de comportement locaux durant leur phase adoucissante.
Les premiers travaux se concentrés sur la mécanique de l’endommagement en quasi
statique puis ont été étendus à la dynamique [MOR 15] pour de la rupture fragile. Des
travaux en parallèle dans [MOË 14] ont permis de coupler l’utilisation d’une vision
locale de l’endommagement dans les zones de la structure n’ayant pas atteintes la phase
adoucissante, et une vision non locale dans les zones de localisation.

Comme on peut le voir figure 1.38, une level set (voir partie [Link]) est utilisée pour
séparer une zone saine d’une zone endommagée. Dans la zone endommagée, la variable
d’endommagement est une fonction explicite de la level set et cette fonction est un pa-
ramètre du modèle. Passée une longueur critique, le matériau est supposé entièrement
endommagé, menant à une transition naturelle vers la rupture. Lorsque l’endommage-
ment atteint 1 dans la zone non locale, la discontinuité des champs cinématiques se met
en place avec la méthode X-FEM, offrant la transition continu/discontinu. C’est ce dernier
point qui différencie cette méthode de l’approche phase field présentée dans le paragraphe
précédent qui n’offre pas cette transition continu/discontinu. Dans [CAZ 15], les auteurs
montrent entre autre que la TLS règle en partie les problèmes de localisation de la fissure
proche d’une surface libre rencontrés avec les phase-fields.
La TLS présente plusieurs avantages par rapport aux méthodes classiques utilisant des
limiteurs de localisation (présentées partie 1.35). D’une part la transition de la localisation
à la rupture est prise en compte dans le modèle. D’autre part, là où les modèles non locaux
en tout point de la structure et à tout instant ne permettent pas de modéliser une dégrada-
tion diffuse, cette méthode permet d’avoir une vision non locale au moment et à l’endroit
nécessaire. Enfin, cela permet de réduire significativement le temps de calcul par rapport
à un modèle entièrement non local. Cette extension a ouvert la voie à la modélisation des
matériaux quasi fragiles et ductiles avec la TLS, notamment en dynamique. Les premiers
travaux à ce sujet sont disponibles dans [MOR 14].

Cette thèse est accessible à l'adresse : [Link]

F IGURE 1.38: Schéma d’une level set épaisse.

1.5 Conclusion
Une liste non exhaustive des différentes modèles et méthodes numériques pour
représenter la propagation de fissure a été présentée. Pour aller plus loin dans la
compréhension des modèles, des ouvrages de référence existent tels que le livre de
Lemaitre et Chaboche [LEM 09] pour les notions d’endommagement dans le contexte
de mécanique des milieux continus, les ouvrages de Bui [BUI 78] ou Anderson
[AND 05] pour les notions de mécanique de la rupture, et enfin l’ouvrage de Freund
[FRE 98] qui s’intéresse plus particulièrement à la mécanique de la rupture en dynamique.

A travers cette section, il ressort qu’il existe principalement deux visions de la

rupture : une vision locale, basée sur la mécanique de l’endommagement, et une vision
énergétique basée sur la mécanique de la rupture. La première permet de modéliser
l’initiation, la croissance et la coalescence de micro fissures pour former la macrofissure.
L’avancée de la fissure se fait "naturellement" grace à la richesse de la loi matériau.
L’autre vision est plus pragmatique et ne décrit pas aussi finement ce qui se passe à
l’échelle locale. Il est alors nécessaire de "rajouter de la physique" via des critères pour
déterminer à la fois vitesse et direction de propagation.

Comme évoqué dans l’introduction générale, les motivations de cette étude sont
de mener des calculs de propagation de fissure en dynamique rapide à l’échelle d’une
structure entière de navire. Or il n’existe pas à ce jour d’outils numériques disponibles
pour simuler une telle propagation.

Nous avons donc fait le choix d’une approche par élément de coque enrichie d’une part
et d’un modèle basé sur la mécanique de la rupture à faible nombre de paramètres d’autre
part, dans l’objectif de maitriser le coût numérique de nos simulation et afin de permettre
leur application à moyen terme sur des structures de grande taille. Ce type d’approche ne
permet pas de simuler la naissance de la fissure et les premiers instant de sa propagation
et nécessite donc l’introduction d’une fissure initiale.

Cette thèse est accessible à l'adresse : [Link]

Vous aimerez peut-être aussi

Analyse des modes de rupture des matériaux
Pas encore d'évaluation
Analyse des modes de rupture des matériaux
46 pages
Mécanique de la Rupture et Fatigue
100% (9)
Mécanique de la Rupture et Fatigue
84 pages
Chap 1
Pas encore d'évaluation
Chap 1
16 pages
Mécanique de la Rupture : Concepts et Modèles
Pas encore d'évaluation
Mécanique de la Rupture : Concepts et Modèles
10 pages
Mécanique de la Rupture : Cours Complet
Pas encore d'évaluation
Mécanique de la Rupture : Cours Complet
54 pages
Loi de Hooke généralisée et rupture des matériaux
Pas encore d'évaluation
Loi de Hooke généralisée et rupture des matériaux
14 pages
Mécanismes de rupture : ductile et fragile
Pas encore d'évaluation
Mécanismes de rupture : ductile et fragile
41 pages
Ruptures ductiles et fragiles en mécanique
Pas encore d'évaluation
Ruptures ductiles et fragiles en mécanique
36 pages
Mécanique de la rupture des matériaux
Pas encore d'évaluation
Mécanique de la rupture des matériaux
20 pages
Ruptures des Liberty Ships en guerre
Pas encore d'évaluation
Ruptures des Liberty Ships en guerre
9 pages
Chap 1
Pas encore d'évaluation
Chap 1
37 pages
Mécanique de la rupture et fatigue des matériaux
Pas encore d'évaluation
Mécanique de la rupture et fatigue des matériaux
16 pages
Meca Rupt Complet New
Pas encore d'évaluation
Meca Rupt Complet New
123 pages
Clément Lemaignan - La Rupture Des Matériaux (2003, EDP Sciences)
100% (1)
Clément Lemaignan - La Rupture Des Matériaux (2003, EDP Sciences)
150 pages
Rupture Des Materiaux
100% (4)
Rupture Des Materiaux
150 pages
Mécanique de la Rupture : Concepts Clés
Pas encore d'évaluation
Mécanique de la Rupture : Concepts Clés
30 pages
Mécanique de la Rupture : Concepts Clés
Pas encore d'évaluation
Mécanique de la Rupture : Concepts Clés
26 pages
Mécanique de la Rupture: Concepts Clés
Pas encore d'évaluation
Mécanique de la Rupture: Concepts Clés
87 pages
Mécanique des fissures dans matériaux élastiques
Pas encore d'évaluation
Mécanique des fissures dans matériaux élastiques
13 pages
Mécanique de Rupture-Abdou-Version 06-12-2023
Pas encore d'évaluation
Mécanique de Rupture-Abdou-Version 06-12-2023
97 pages
Mécanique de la rupture et fatigue
Pas encore d'évaluation
Mécanique de la rupture et fatigue
10 pages
Mécanique de la Rupture et Matériaux
Pas encore d'évaluation
Mécanique de la Rupture et Matériaux
19 pages
Faciès de Rupture des Matériaux
Pas encore d'évaluation
Faciès de Rupture des Matériaux
4 pages
Deux PDF
Pas encore d'évaluation
Deux PDF
15 pages
Mécanique de la Rupture et Fatigue
100% (1)
Mécanique de la Rupture et Fatigue
46 pages
Mécanique de la rupture : Historique et concepts
88% (8)
Mécanique de la rupture : Historique et concepts
17 pages
Fiabilité des matériaux et rupture mécanique
Pas encore d'évaluation
Fiabilité des matériaux et rupture mécanique
33 pages
Chap 1
Pas encore d'évaluation
Chap 1
24 pages
Historique de la Mécanique de Rupture
Pas encore d'évaluation
Historique de la Mécanique de Rupture
10 pages
Analyse des faciès de rupture
Pas encore d'évaluation
Analyse des faciès de rupture
3 pages
Rupture des matériaux : types et mécanismes
Pas encore d'évaluation
Rupture des matériaux : types et mécanismes
22 pages
Introduction à la Mécanique de la Rupture
100% (2)
Introduction à la Mécanique de la Rupture
47 pages
Rupture mécanique : ductile et fragile
100% (1)
Rupture mécanique : ductile et fragile
7 pages
Analyse de la rupture du bitume par FEM
Pas encore d'évaluation
Analyse de la rupture du bitume par FEM
7 pages
Fatigue des structures aérodynamiques
Pas encore d'évaluation
Fatigue des structures aérodynamiques
15 pages
Mécanique de la Rupture : Concepts Clés
100% (2)
Mécanique de la Rupture : Concepts Clés
20 pages
Chapitre 1 MR
Pas encore d'évaluation
Chapitre 1 MR
18 pages
Mécanique de la Rupture : Concepts Clés
100% (5)
Mécanique de la Rupture : Concepts Clés
156 pages
Mécanique de la rupture et essais associés
100% (5)
Mécanique de la rupture et essais associés
35 pages
Mécanique de la rupture par fatigue
Pas encore d'évaluation
Mécanique de la rupture par fatigue
127 pages
Rupture Fragile Des Matériaux SM - YANE
Pas encore d'évaluation
Rupture Fragile Des Matériaux SM - YANE
30 pages
Devis pour machines d'injection RTM
Pas encore d'évaluation
Devis pour machines d'injection RTM
75 pages
Introduction
Pas encore d'évaluation
Introduction
1 page
Mécanique de la Rupture et Endommagement
Pas encore d'évaluation
Mécanique de la Rupture et Endommagement
42 pages
Analyse du facteur d'intensité de contrainte
Pas encore d'évaluation
Analyse du facteur d'intensité de contrainte
63 pages
Introduction
Pas encore d'évaluation
Introduction
20 pages
Mécanique de la rupture des matériaux
Pas encore d'évaluation
Mécanique de la rupture des matériaux
173 pages
Mécanique de la Rupture : Concepts Clés
Pas encore d'évaluation
Mécanique de la Rupture : Concepts Clés
32 pages
Mécanique de la Fracture Élastique
Pas encore d'évaluation
Mécanique de la Fracture Élastique
90 pages
Analyse des faciès de rupture des matériaux
Pas encore d'évaluation
Analyse des faciès de rupture des matériaux
11 pages
Admin,+11+-Pluvinage C4
Pas encore d'évaluation
Admin,+11+-Pluvinage C4
5 pages
Types et Mécanismes de Rupture des Matériaux
Pas encore d'évaluation
Types et Mécanismes de Rupture des Matériaux
45 pages
Poly-Rup - PDF 5
Pas encore d'évaluation
Poly-Rup - PDF 5
137 pages
Déformation et rupture des matériaux
Pas encore d'évaluation
Déformation et rupture des matériaux
9 pages
Mécanique de la Rupture des Matériaux
100% (1)
Mécanique de la Rupture des Matériaux
13 pages
Mécanique de la rupture : types et facteurs
Pas encore d'évaluation
Mécanique de la rupture : types et facteurs
15 pages
Mécanique de la Rupture : Concepts Clés
Pas encore d'évaluation
Mécanique de la Rupture : Concepts Clés
12 pages
Mécanique de La Rupture: Jean-Luc ENGERAND
Pas encore d'évaluation
Mécanique de La Rupture: Jean-Luc ENGERAND
12 pages
Mécanique de La Rupture
100% (2)
Mécanique de La Rupture
12 pages
Guide méthodologique d'étude de dangers
100% (2)
Guide méthodologique d'étude de dangers
19 pages
Pénsum Génie Civil : Cours et Crédits
Pas encore d'évaluation
Pénsum Génie Civil : Cours et Crédits
1 page
Citations marquantes d'Albert Einstein
Pas encore d'évaluation
Citations marquantes d'Albert Einstein
3 pages
Dimensionnement d'un système de pompage solaire
Pas encore d'évaluation
Dimensionnement d'un système de pompage solaire
22 pages
Intégrales et Méthodes de Calcul en TS2A
Pas encore d'évaluation
Intégrales et Méthodes de Calcul en TS2A
2 pages
Mécanique du Solide : Exercices et Corrections
Pas encore d'évaluation
Mécanique du Solide : Exercices et Corrections
2 pages
Calculs de RDM pour tube rond entre appuis
Pas encore d'évaluation
Calculs de RDM pour tube rond entre appuis
4 pages
Modélisation d'un moteur asynchrone 6 phases
Pas encore d'évaluation
Modélisation d'un moteur asynchrone 6 phases
4 pages
Devoir de synthèse 3ème Eco-G 2024-2025
Pas encore d'évaluation
Devoir de synthèse 3ème Eco-G 2024-2025
2 pages
Optimisation Déterministe avec Contraintes
Pas encore d'évaluation
Optimisation Déterministe avec Contraintes
19 pages
Diagrammes de Rayonnement Antennes MATLAB
Pas encore d'évaluation
Diagrammes de Rayonnement Antennes MATLAB
2 pages
Analyse des données MIDO 2020-2021
Pas encore d'évaluation
Analyse des données MIDO 2020-2021
21 pages
Traction électrique et locomotives SNCF
Pas encore d'évaluation
Traction électrique et locomotives SNCF
38 pages
Feuille d'exercices sur les nombres complexes
Pas encore d'évaluation
Feuille d'exercices sur les nombres complexes
3 pages
Planning des examens 2023-2024 Maroc
Pas encore d'évaluation
Planning des examens 2023-2024 Maroc
1 page
Lois de pilotage en mécanique rotative
Pas encore d'évaluation
Lois de pilotage en mécanique rotative
30 pages
Volet Bois Cadre Lames à la Française 34 mm
Pas encore d'évaluation
Volet Bois Cadre Lames à la Française 34 mm
90 pages
Normes des Thrombocytes en Physique
Pas encore d'évaluation
Normes des Thrombocytes en Physique
10 pages
Mesure du Champ Magnétique par Effet Hall
100% (1)
Mesure du Champ Magnétique par Effet Hall
5 pages
Liste de pièces pour Junkers ZWE 24-4
Pas encore d'évaluation
Liste de pièces pour Junkers ZWE 24-4
22 pages
Trajectoire d'une boule de pétanque
Pas encore d'évaluation
Trajectoire d'une boule de pétanque
6 pages
Agronomie I : Eau - Support de cours
Pas encore d'évaluation
Agronomie I : Eau - Support de cours
16 pages
Symétrie des fonctions et courbes
Pas encore d'évaluation
Symétrie des fonctions et courbes
2 pages
Énergie cinétique et potentielle en physique
Pas encore d'évaluation
Énergie cinétique et potentielle en physique
2 pages
Rapport de Stage chez Tunicast
100% (1)
Rapport de Stage chez Tunicast
41 pages
TD Circuits Magnétiques : Exercices Corrigés
Pas encore d'évaluation
TD Circuits Magnétiques : Exercices Corrigés
3 pages
Exercices corrigés d'algèbre 2
Pas encore d'évaluation
Exercices corrigés d'algèbre 2
2 pages
Introduction à l'Électromagnétisme
Pas encore d'évaluation
Introduction à l'Électromagnétisme
147 pages
Introduction à la convection en MEC3200
Pas encore d'évaluation
Introduction à la convection en MEC3200
36 pages