Chapitre 2

Extension de l’approche X-FEM en

dynamique rapide pour la propagation
de fissure dans des coques de type
Mindlin-Reissner

Sommaire
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.2 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner 52
2.2.1 Objectifs et hypothèses . . . . . . . . . . . . . . . . . . . . . . . . 52
2.2.2 Propriétés de l’élément fini Q4GS . . . . . . . . . . . . . . . . . . 53
2.2.3 Schéma temporel d’intégration . . . . . . . . . . . . . . . . . . . . 57
2.2.4 Création de l’élément de coque enrichi X-Q4GS. . . . . . . . . . . 61
2.2.5 Description de la fissure . . . . . . . . . . . . . . . . . . . . . . . 66
2.2.6 Cas test à fissure fixe . . . . . . . . . . . . . . . . . . . . . . . . . 70
2.3 Critère de propagation de fissures en dynamique explicite . . . . . . . . 74
2.3.1 Choix d’un critère de propagation . . . . . . . . . . . . . . . . . . 74
2.3.2 Direction de propagation en traction/cisaillement . . . . . . . . . . 76
2.3.3 Incrément de propagation dans un pas de temps . . . . . . . . . . . 83
2.3.4 Extension du critère pour les coques . . . . . . . . . . . . . . . . . 85
2.3.5 Moyennes temporelles sur les avancées du front de fissure . . . . . 87
2.3.6 Mise en œuvre numérique de la propagation . . . . . . . . . . . . . 89
2.4 Simulation d’un exemple type issu de la littérature . . . . . . . . . . . . 97
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

51

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

2.1 Introduction
La méthode des éléments finis étendus développée au début des années 2000
[MOË 99] commence à avoir une certaine maturité et apparaît peu à peu dans les codes
de calcul utilisés dans l’industrie. Ce modèle basé sur la partition de l’unité [MEL 96]
consiste à enrichir la base polynomiale d’approximation de manière à rendre possible
la prise en compte entre autres d’une discontinuité au sein d’un élément. De nombreux
travaux ont été menés pour traiter des cas de propagation de fissure en dynamique
rapide sur des modèles 2D plans ou 3D volumiques. Cependant, le maillage doit être
relativement fin pour décrire correctement les différents champs dans le voisinage de la
fissure pour ensuite être utilisés dans un critère de propagation. Mener de tels calculs à
l’échelle d’une structure navale dans un contexte de bureau d’étude paraît délicat. Ces
structures étant élancées (assemblages de plaques et de coques), il est pertinent d’utiliser
des modèles coques qui permettent de diminuer considérablement la taille d’un modèle
et de gagner en temps de calcul.

Ce chapitre propose de détailler les différentes étapes nécessaires afin de pouvoir

simuler la propagation de fissure en dynamique rapide dans des éléments finis de coques
à l’aide de la méthode des éléments finis étendus.

La première partie de ce chapitre est consacrée à la description de tous les "ingré-

dients" qu’il est nécessaire d’ajouter à l’élément fini de coque afin que les champs d’in-
connues cinématiques de celui-ci puissent décrire spatialement la fissure. La seconde se
concentre davantage sur le critère de propagation à savoir la détermination des grandeurs
permettant de définir la propagation, en terme de vitesse et de direction.

2.2 Couplage de la méthode X-FEM au modèle coque de

Mindlin-Reissner
2.2.1 Objectifs et hypothèses
Certains travaux traitent le cas de plaques fissurées en flexion [LAS 09] ou de coques
fissurées [ARE 05] en se basant sur la théorie de Kirchhoff-Love des coques fines.
Toutefois, les épaisseurs conséquentes des coques qui forment des navires comme les
sous-marins imposent le choix d’une autre formulation plus adaptée aux coques dites
"épaisses" ; c’est la théorie de Mindlin-Reissner. Celle-ci permet la prise en compte du
cisaillement transverse mais nécessite un traitement particulier pour éviter le phénomène
dit de "verrouillage" lorsque l’élancement (longueur/épaisseur) de la coque devient
trop important. La méthode X-FEM a été utilisée avec cette théorie par Dolbow et al.
[DOL 00b] à l’aide de l’élément MITC4 pour traiter le cas des plaques en élasticité
bidimensionnelle.

52

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

Afin de traiter les cas de coque tridimensionnelle fissurée basée sur la théorie de
Mindlin-Reissner, le choix a été fait d’enrichir le quadrangle Q4γ24. Cet élément est
brièvement présenté dans le paragraphe 2.2.2

Des essais de détonation ont été effectués sur des éprouvettes (plates et cylindriques)
pré-entaillées sur le site de Tourris, en collaboration avec DGA Techniques Navales.
Lorsque le défaut initial n’est pas débouchant, les essais montrent que la fissure se propage
dans un premier temps dans l’épaisseur de l’éprouvette. La fissure se propage ensuite en
restant traversante. Lors de la phase de propagation, la fissure traversante se propage avec
une inclinaison d’environ 45 degrés par rapport au plan moyen de la coque ce qui traduit
une rupture en présence d’un cisaillement. De plus, le phénomène de striction apparait en
pointe de fissure, qui se traduit par une réduction de l’épaisseur de la coque localement.
On constate une inclinaison de 45 degrés lorsqu’il s’agit de cisaillement pur et d’un angle
légèrement inférieur pour un chargement mixte comme on peut le remarquer sur la figure
2.1. A noter que lorsqu’une entaille introduite est initialement débouchantes, le même
type de propagation est observé.

F IGURE 2.1: Orientation de la fissure dans l’épaisseur de la structure.

Les modèles de coques étant formulés à partir de la surface moyenne de la structure,

l’analyse de l’évolution de la fissure à travers l’épaisseur n’est guère prévisible. De plus,
comme il l’a été rappelé en partie [Link], la méthode X-FEM ne permet pas de traiter la
phase d’amorçage de fissures. Ce travail se focalise donc sur des fissures préexistantes qui
seront supposées débouchantes dans l’épaisseur.

2.2.2 Propriétés de l’élément fini Q4GS

Le Q4GS est un élément implémenté dans Europlexus qui permet d’obtenir les so-
lutions numériques des problèmes de plaque et coque, qu’elles soit minces ou épaisses.
Le cisaillement transverse est pris en compte avec les hypothèses de Mindlin-Reissner.
Pour cet élément, la normale en tout point est obtenue par interpolation des normales aux

53

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

nœuds. Il s’agit de la vision coque dégénérée définie dans la partie [Link]. C’est cette
version simplifiée de coque épaisse que l’on propose d’enrichir avec la méthode X-FEM
dans la suite de nos travaux.

(-1,1) (1,1)
4 3

1 2
(-1,-1) (1,-1)

F IGURE 2.2: Elément de référence du quadrangle à 4 nœuds.

C’est un quadrangle à épaisseur constante dont la partie membrane est représentée par
l’élément Q4 et la partie flexion/cisaillement par l’élément Q4γ (F IG .2.2).

[Link] Caractéristiques
C’est un élément courbe à quatre nœuds et six degrés de liberté par nœud. Ses incon-
nues nodales sont les déplacements locaux ux , uy et uz , et les rotations locales θx , θy et
θz ;
– ux et uy représentent les effets de membrane,
– θx et θy représentent les effets de flexion,
– uz , θx et θy servent à établir les effets de cisaillement transversal,
– θz est utilisée pour éviter la singularité de la matrice de rigidité.

Les variables nodales après assemblage sont les déplacements UX , UY et UZ des trois
noeuds et les rotations θX , θY et θZ autour des axes globaux X, Y et Z (voir F IG .2.3).
L’intégration numérique se fait sur quatre points dans le plan de l’élément de référence
(voir TAB .2.1), et cinq points dans l’épaisseur en faisant varier z entre − 2h et h2 .

Alors que les trois nœuds d’un triangle sont nécessairement coplanaires, ce n’est plus
le cas pour les quatre nœuds d’un quadrangle. Il est alors nécessaire de définir un repère
local en chaque sommet et en chaque point d’intégrations pour prendre en compte la

54

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

F IGURE 2.3: Variables nodales locales (5 ddl/noeud).

courbure de l’élément.

[Link] Fonctions d’interpolation

Que ce soit pour les déplacements ou les rotations, les fonctions de forme classiques
bilinéaires des quadrangles sont utilisées :

1
N1 (ξ, η) = ((1 − ξ)(1 − η)) (2.1)
4
1
N2 (ξ, η) = ((1 + ξ)(1 − η)) (2.2)
4
1
N3 (ξ, η) = ((1 + ξ)(1 + η)) (2.3)
4
1
N4 (ξ, η) = ((1 − ξ)(1 + η)) (2.4)
4

[Link] Traitement du verrouillage en cisaillement transverse

La prise en compte du cisaillement transverse nécessite quelques précautions car il
entraine quelques complications avec notamment l’apparition du verrouillage. Pour le
Q4GS, c’est la méthode du champ assumé de déformation de cisaillement transverse

55

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

TABLEAU 2.1: Points d’integration du Q4GS.

Coordonnées Poids
ξi ηi ωi

√ √
1/ 3 1/ 3
√ √
−1/ 3 1/ 3 1
√ √
−1/ 3 −1/ 3
√ √
1/ 3 −1/ 3

qui est utilisée (F IG .2.4). Le cisaillement transversal γ n’est pas directement calculé aux
quatre points d’intégration comme pour la membrane et la flexion, mais il est interpolé
à l’aide des valeurs calculées aux quatre points (A1, A2, B1 et B2) milieux des côtés du
quadrilatère. Ainsi, avec [J] la matrice jacobienne de la transformation, le cisaillement
transverse s’interpole de la manière suivante :
 1−η
A1 + 1+η γA2 

      2 ξ
 γ 2 ξ  
γxz −1 γξz
= [J] = (2.5)
γyz γηz  
 1−ξ B1 1+ξ B2 
 
2 γη + 2 γη
avec γξz et γηz constants sur chaque côté de l’élément. Ce traitement du cisaillement
permet d’éviter le phénomène de verrouillage en cisaillement transverse.

4
A2 3

B1 B2

A1
2
1

F IGURE 2.4: Champ assumé de déformation en cisaillement transverse.

56

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

Une bonne description de l’élément peut être trouvée dans Batoz [BAT 92] et Zeng
[ZEN 98].

2.2.3 Schéma temporel d’intégration

Dans cette partie on s’intéresse à la manière d’intégrer numériquement l’équation
d’équilibre dynamique. On présente en premier le cas général d’un problème sans amor-
tissement, puis les modifications engendrée par l’introduction d’un amortissement.

[Link] Sans amortissement

Pour le cas d’un problème sans contact ni amortissement, l’équation d’équilibre dis-
crétisée peut s’écrire de la manière suivante :

MÜ = Fext − Fint , (2.6)

où Fext sont les efforts extérieurs et Fint les efforts intérieurs. Dans le cas élastique linéaire
on a laRrelation Fint = KU et dans le cas non linéaire matériau on utilise généralement
Fint = Ω BT σdΩ où Ω est la configuration déformée, σ la contrainte de Cauchy et B
l’opérateur divergence discrétisé.

Deux grands types de schéma temporels d’intégration existent pour cela : les schémas
explicites et implicites en temps. Les premiers permettent une résolution rapide mais
sont conditionnellement stable. Cela implique un pas de temps suffisamment petit pour
éviter que le calcul ne diverge. Les schémas d’intégration implicites permettent d’utiliser
un pas de temps qui peut être important mais nécessitent des itérations pour assurer la
convergence de l’équilibre à chaque pas de temps, rendant chaque itération temporelle
plus longue à résoudre.

En dynamique rapide, les problèmes à traiter sont des cas de chargements brefs et il
est en général préférable d’avoir des pas de temps suffisamment petits afin de suivre au
mieux les variations rapides du chargement. Ainsi, le schéma implicite est souvent écarté
car la durée des pas de temps imposée par la physique du problème conduiraient à de trop
grands temps de calcul.

L’algorithme d’intégration en temps utilisé dans Europlexus est de type Newmark

[NEW 59]. Il est basé sur un développement limité de Taylor des fonctions temporelles.
Les dérivations successives peuvent être approchées par un développement de Taylor en
utilisant deux paramètres, β pour le déplacement U et γ pour la vitesse U̇. En notant Un
le déplacement à l’instant n, U̇n la vitesse et Ün l’accélération, le passage du pas n au pas

57

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

n + 1 est donné par :

∆t 2 ∆t 3
 
Ün+1 − Ün
Un+1 = Un + ∆tU̇n + Ün + 3β , (2.7)
2 3 ∆t
Ün+1 = M −1 (Fext,n+1 − Fint,n+1 ), (2.8)
∆t 2
 
Ün+1 − Ün
U̇n+1 = U̇n + ∆tÜn + 2γ . (2.9)
2 ∆t
En réarrangeant les équations, il vient :
1
Un+1 = Un + ∆tU̇n + ∆t 2 ( − β)Ün + ∆t 2 βÜn+1 , (2.10)
2
−1
Ün+1 = M (Fext,n+1 − Fint,n+1 ), (2.11)
U̇n+1 = U̇n + ∆t(1 − γ)Ün + γ∆tÜn+1 . (2.12)

Le tableau TAB .2.2 présente les différents algorithmes utilisés en fonction du couple
de paramètres (β, γ).

TABLEAU 2.2: Résumé des différents schémas numériques de Newmark.

Algorithme Type γ β
Purement explicite Explicite 0 0
Différence centrée Explicite 1/2 0
Accélération linéaire Implicite 1/2 1/6
Accélération moyenne Implicite 1/2 1/4
Accélération moyenne modifiée Implicite 1/2 + γ (1 + γ)2 /4
Fox-Goodwin Implicite 1/2 1/2

Le schéma utilisé dans Europlexus est donné pour le couple de valeur (β, γ) = (0, 21 ).
Il s’agit du schéma explicite de différence centrée, dont les relations simplifiées sont don-
nées par :
∆t 2
Un+1 = Un + ∆tU̇n + Ün , (2.13)
2
Ün+1 = M −1 (Fext,n+1 − Fint,n+1 ), (2.14)
∆t ∆t
U̇n+1 = U̇n + Ün + Ün+1 . (2.15)
2 2
Ainsi l’algorithme permet une grande souplesse pour traiter les problèmes de
dynamique rapide et les comportements non linéaires. Le déplacement au pas de temps
tn+1 ne dépend que des éléments de solution au temps tn ce qui rend la résolution directe.

58

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

[Link] Avec amortissement

Lorsqu’il y a de l’amortissement ou du contact, un terme dépendant de la vitesse U̇ ap-
paraît dans l’équation d’équilibre 2.6. Si on note C la matrice d’amortissement, l’équation
devient :
MÜ +CU̇ = Fext − Fint . (2.16)
Le passage du pas n au pas n + 1 est alors donné par :
∆t 2
Un+1 = Ün ,
Un + ∆tU̇n + (2.17)
2
Ün+1 = M −1 (Fext,n+1 − Fint,n+1 −CU̇n+1 ), (2.18)
∆t ∆t
U̇n+1 = U̇n + Ün + Ün+1 . (2.19)
2 2
On remarque alors que le schéma est explicite en déplacement, mais qu’il ne l’est
plus en vitesse et en accélération : l’accélération au temps n + 1 dépend désormais de
la vitesse au temps n + 1 et inversement. Ainsi la résolution n’est plus directe et des
itérations sont nécessaires.

Pour que l’algorithme de la différence centrée reste

explicite, Belytschko [BEL 99]
propose d’introduire la vitesse au pas de temps n + 12 définie de la manière suivante :
Un+1 −Un
U̇n+ 1 =. (2.20)
2 ∆t
En injectant l’équation 2.17 dans l’équation 2.20 il vient :
∆t
U̇n+ 1 = U̇n + Ün . (2.21)
2 2
Les relations 2.13 et 2.15 écrites au pas de temps n, la définition de U̇n− 1 et l’équation
2
2.21 permettent d’écrire :

U̇n+ 1 = U̇n− 1 + ∆tÜn . (2.22)

2 2

L’introduction de U̇n+ 1 dans le schéma de la différence centrée modifie l’équation

2
d’équilibre qui s’écrit alors à l’instant n+1 :
MÜn+1 +CU̇n+ 1 = Fext,n+1 − Fint,n+1 . (2.23)
2

Finalement, les équations correspondantes à ce schéma d’intégration s’écrivent :

∆t
U̇n+ 1 = U̇n + Ün , (2.24)
2 2
∆t 2
Un+1 = Un + ∆tU̇n + Ün , (2.25)
2
Ün+1 = M −1 (Fext,n+1 − Fint,n+1 −CU̇n+ 1 ), (2.26)
2
∆t ∆t
U̇n+1 = U̇n + Ün + Ün+1 , (2.27)
2 2

59

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

ce qui permet d’obtenir à nouveau un schéma explicite. L’initialisation de la vitesse se fait

alors de la manière suivante :
∆t
U̇ 1 = U̇0 + Ü0 . (2.28)
2 2

[Link] Pas de temps critique pour les schémas explicites

La condition de stabilité d’un schéma explicite pour le cas unidimensionnel s’exprime

de la manière suivante :
∆x
∆t 6 ∆tcrit = , (2.29)
c

avec ∆x la plus petite longueur

q caractéristique de la structure discrétisée, et c la célérité
de l’onde donnée par c = Eρ où E est le module d’Young et ρ la masse volumique.
Physiquement, ∆tcrit correspond au temps mis par l’onde de célérité c pour traverser
la plus petite longueur de l’élément le plus petit. Sur l’exemple d’une structure unidi-
mensionnelle figure 2.5, on peut comprendre ce critère en imaginant que l’information se
transmet de la gauche vers la droite :
– Si ∆t 6 ∆tcrit les phénomènes peuvent être décris dans tous les éléments,
– Si ∆t > ∆tcrit certains éléments ne voient pas l’information.

F IGURE 2.5: Stabilité des calculs en fonction du pas de temps.

Il y a alors stabilité si l’incrément de temps permet de décrire l’information dans tous

les éléments constituant la structure. Il faut donc veiller à choisir subtilement la taille
de ses éléments car quelques petits éléments pour capter un gradient de contrainte vont
pénaliser le temps de calcul pour toute la structure. Un schéma implicite de Newmark est
lui toujours stable quelque soit le pas de temps tant que 2β ≥ γ ≥ 21 .

60

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

2.2.4 Création de l’élément de coque enrichi X-Q4GS.

[Link] Enrichissement de l’élément

Concernant l’enrichissement, seules les fonctions sauts de Heaviside sont utilisées

pour décrire la discontinuité à la fois en déplacement (2.30) et en rotation (2.31). Il a
été montré dans [GRE 07] que les fonctions singulières en pointe de fissure ne sont pas
indispensables pour mieux décrire le champ de déplacement, pourvu que le maillage soit
suffisamment fin au voisinage du front.

ui (x) = ∑ NA (x)uiA + ∑cut NB(x)H (x)u∗iB, (2.30)

A∈N B∈N

θi (x) = ∑ NA (x)θiA + ∑cut NB(x)H (x)θ∗iB, (2.31)

A∈N B∈N


1 si x est au dessus de la fissure
avec H (x) = (2.32)
−1 si x est en dessous de la fissure

Il convient de remarquer qu’avec ce choix d’enrichissement, la position de la pointe

n’est qu’approximatif étant donné que la propagation a lieu sur un nombre entier d’élé-
ments. Dès que le calcul prévoit une avancée de la pointe de fissure dans un nouvel élé-
ment, celui-ci est entièrement tranché comme on peut le voir figure 2.8.

F IGURE 2.6: Stratégie d’enrichissement pour une fissure quelconque placée sur un
maillage.

61

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

Diagonalisation de la matrice masse Pour profiter pleinement du schéma temporel

explicite et ne travailler que sur une résolution vectorielle, la matrice masse doit être
diagonalisée. Ainsi, l’équation d’équilibre peut s’exprimer sous la forme :

MdiagÜ = Fext − Fint , (2.33)

De base, la matrice de masse consistante prend la forme générale suivante avec un

rangement adéquat des fonctions de forme :
 
Mstandard Mcouplage
MXFEM = , (2.34)
Mcouplage Menrichie

où la matrice Mstandard est la matrice de masse construite à partir des fonctions de forme
non enrichies seulement, la matrice Menrichie construite à partir des fonctions de forme
enrichies uniquement, et Mcouplage par le couplage des fonctions de forme standard et
enrichies.
La méthode de diagonalisation de matrice masse basée sur la conservation de
l’énergie cinétique mise en place par Menouillard et al. [MEN 07] a été reprise. Les
auteurs montrent en outre que les masses associées aux degrés de liberté enrichis peuvent
prises égales à celles des degrés de liberté classiques, tout en assurant la conservation de
l’énergie cinétique pour des mouvements de corps rigide qu’il y ait une fissure ou non.

Dans le cas des coques X-FEM, on choisit alors non seulement de prendre égales
les masses associées aux degrés de liberté classiques et enrichis en déplacement, mais
également pour les masses associées aux degrés de liberté classiques et enrichis relatives
aux rotations. La matrice masse diagonale de taille (12 x12) d’un élément de coque enrichi
devient donc :
!
diag
diag Mstandard 0
MXFEM = diag , (2.35)
0 Menrichie
avec : !
f em
diag mu 0
Mstandard = f em , (2.36)
0 mθ
!
x f em
diag mu 0
Menrichie = x f em , (2.37)
0 mθ

en prenant :
x f em
muf em = mu∗ , (2.38)
f em x f em
mθ = mθ∗ . (2.39)

62

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

Pas de temps critique de l’élément enrichi Gerlach [GER 99] a montré que l’en-
richissement discontinu introduit un caractère non stabilisant du schéma de Newmark
explicite. Comme expliqué précédemment, le pas de temps critique dépend de la plus
petit taille d’élément. Ainsi pour un élément dont un des côté est coupé, c’est la plus
petite des parties du côté coupé qu’il faut prendre en compte pour le calcul du pas de
temps, si bien que le pas de temps explicite tend vers zéro lorsque la fissure coupe un
élément très près d’un noeud.

Menouillard et al. ont montré dans [MEN 06] puis [MEN 08] que pour le cas d’ élé-
ments poutres ou triangles, l’utilisation de la matrice de masse diagonalisée, souvent ap-
pelée matrice lumpée, conduisait à un pas de temps critique non nul. Ils ont même prouvé
qu’en utilisant une matrice diagonale par bloc le pas de temps critique est le même que si
f em
l’élément n’était pas coupé. Ces résultats sont regroupés dans TAB .2.3 en appelant ∆tc
le pas de temps critique de l’élément non coupé pour lequel on utilise la matrice de masse
x f em
diagonale, et ∆tc le pas de temps critique de l’élément coupé.
TABLEAU 2.3: Valeur minimum du pas de temps critique normalisé en fonction de
la matrice masse utilisée pour un élément X-FEM - Éléments poutres et triangles.

x f em f em
Matrice masse ∆tc /∆tc
Consistente 0√
Diagonale 1/ 2
Diagonale par bloc 1

Le cas des quadrangles et des éléments volumiques sont traités dans [MEN 07,
MEN 08],. Les auteurs montrent qu’avec toutes les méthodes de diagonalisation et quel
que soit le type d’élément fini étudié, un pas de temps de √12 fois le pas de temps de l’élé-
ment sans enrichissement garantit la stabilité. En pratique, le pas de temps critique choisi
pour assurer la stabilité du calcul est est donné par :
f em
∆tc
∆tcx f em = . (2.40)
2

[Link] Technique d’intégration

Numériquement, la prise en compte de l’enrichissement intervient dans le calcul des
forces intérieures de l’élément fini. Il faut en effet venir ajouter les contributions des
degrés de liberté enrichis aux noeuds concernés par un enrichissement actif. L’expression
des forces intérieures devient alors :
Z Z
Fint = BT .σdΩ + HBT .σdΩ . (2.41)
| Ω {z } | Ω {z }
Fstandard Fenrichie

63

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

La présence de fonctions discontinues au sein des éléments impose de modifier les

règles d’intégration numérique. Le champ de déplacement n’est en effet plus polynomial
du fait de la fonction saut, ce qui ne permet plus d’utiliser une quadrature classique.

Une première approche développée par Moës et al. [MOË 99] est de sous découper
les éléments tranchés de manière compatible avec la fissure, et ainsi intégrer exactement
des fonctions continues de part et d’autre de la fissure. Ces sous-découpages en triangles
n’interviennent dans le calcul que pour l’intégration numérique. Il ne s’agit en aucun cas
de nouveaux éléments plus fins mais de géométries sur lesquelles les règles d’intégrations
classiques sont connues.
Cette méthode parfaitement adaptée aux matériaux élastiques peut néanmoins
entrainer une non conservation de l’énergie pour des matériaux non linéaires. En effet,
l’utilisation d’un telle technique conduit à un changement du nombre et de la position
des points de Gauss lorsque la fissure se propage. Il faut dans ce cas effectuer une
projection des champs élémentaires lors de l’étape de remaillage des sous-éléments
(contraintes, déformations etc.), dans une zone où leur gradient est particulièrement élevé
(voir F IG .2.7).

(a) Avant propagation. (b) Après propagation.

F IGURE 2.7: Intégration par sous découpage, projection des champs sur de
nouveaux points de Gauss.

Pour éviter ces problèmes et l’étape de sous-découpage des éléments, Elguedj

[ELG 06b] proposent de subdiviser au préalable par un nombre fixe de sous-éléments
de forme similaire tous les éléments enrichis. La position des points d’intégrations est
par conséquent connue avant même que la plasticité (entre autres) ne se développe dans
l’élément. Il n’y a alors plus d’étape de projection. L’intégration n’est plus exacte mais
l’erreur est minimisée en utilisant un nombre accru de sous éléments et les contributions
sont pondérées par l’aire de part et d’autre de la fissure. Pelée de Saint Maurice [SAI 14]
a montré que pour minimiser l’erreur, l’important n’est pas tant le nombre de sous
éléments mais plutôt le nombre de points d’intégration qui interviennent dans chaque
zone. La figure 2.8 présente le cas d’un sous-découpage en 4×4 sous éléments ayant

64

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

chacun quatre points d’intégration dans le plan, soit un total de 64 points de Gauss.

F IGURE 2.8: Intégration numérique par sous-découpage régulier

Une autre solution intéressante pour traiter l’intégration numérique des matériaux
non linéaires a été proposée par Ventura [VEN 06]. Dans les éléments coupés par la
fissure, l’enrichissement par la fonction de Heaviside est remplacé par une nouvelle
fonction continue polynomiale dont le degré dépend du type d’élément et de la position
de l’élément. Cela permet d’une part d’éviter l’étape de projection des champs et
également de réduire le nombre de points d’intégration par rapport à la méthode dé-
crite précédemment. Récemment, Martin et al. [MAR 15] ont automatisé le calcul de la
fonction d’enrichissement polynomiale pour des éléments bi et tridimensionnels linéaires.

La simplicité d’implémentation et les bons résultats affichés par la méthode de

découpage en sous-éléments proposée par Elguedj [ELG 06b] (notamment pour des
matériaux élastoplastiques) nous ont conduit à choisir cette méthode d’intégration.

Remarque : Dans ses travaux, Elguedj utilise la méthode d’intégration à l’aide de sous
éléments uniquement dans une zone située dans un rayon R plastique centrée sur la pointe
de fissure. Lorsque la plasticité est confinée, cela permet en effet de réduire le nombre
d’intégration du problème. Dans les essais dynamiques que nous serons amenés à si-
muler, la plasticité de se limite pas uniquement à la pointe de fissure et c’est pourquoi
cette stratégie n’est pas mise en place. Ainsi l’ensemble des éléments X-FEM du maillage

65

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

utilise une intégration à l’aide de sous-éléments.

2.2.5 Description de la fissure

[Link] Choix de la description
L’expérience montre qu’en général une fissure présente sur une des peaux de la coque
se propage dans un premier temps dans l’épaisseur pour devenir débouchante, puis se
propage à l’ensemble de la coque. C’est pourquoi on rappelle que l’hypothèse de fissure
traversante dans l’épaisseur a été faite.

Une méthode couramment employée pour pour représenter la fissure est l’utilisation
de fonctions de niveau appelées "level-sets" en anglais, proposée par Sethian [SET 96],
puis couplé à la méthode X-FEM par Stolarska et al. [STO 01]. Ces deux fonctions de
niveau (φ,ψ) définissent implicitement la position du plan et du front de fissure comme
on peut le voir figure 2.9. Gravouil et al. [GRA 02] ont ensuite généralisé l’utilisation des
fonctions de niveau au cas tridimensionnel.

F IGURE 2.9: Représentation de la fissure sous forme de segments sur le plan

moyen de la coque

Mathématiquement, la fissure est représentée par l’ensemble des points tels que φ = 0
et ψ 6 0. Le front de fissure est alors décrit par φ = 0 et ψ = 0. L’avantage principal
de cette description est d’avoir accès à la distance et la position relative de n’importe
quel point par rapport à la fissure. Cela est notamment intéressant pour affecter le signe
de l’enrichissement saut. Cependant, l’étape d’actualisation de ces level-sets par la
résolution des équations d’Hamilton-Jacobi ou par méthode géométrique, ainsi que leurs
ré-initialisation et ré-orthogonalisation posent aujourd’hui encore quelques problèmes de

66

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

robustesse, surtout dans le cas de grands déplacements et de grandes distorsions comme

l’a montré Pelée de Saint Maurice [SAI 14].

La description choisie ici se base sur le fait que la fissure est contenue dans le
plan moyen de la coque. Ainsi, la géométrie de la coque donne déjà deux informations
spatiales pour localiser la fissure. L’utilisation d’un maillage d’éléments unidimension-
nels sur le plan moyen de la coque présentée figure 2.10 permet donc de connaitre
complètement la position de la fissure. Avec cette méthode, la fissure est représentée
explicitement. Son principal avantage est la facilité de mettre à jour la fissure. Néanmoins
la gestion de la normale au front est moins naturelle qu’avec les level-sets et l’affectation
des enrichissement n’est pas aussi direct.

Ce type de méthode consistant à utiliser un maillage pour représenter explicitement

la fissure connait un regain d’intérêt ces dernières années avec notamment l’approche
implicites/explicite proposée par Fries et al. [FRI 12] ou Prabel et al. [PRA 11]. L’idée est
de conserver la level set pour renseigner où et comment enrichir et d’utiliser un maillage
(uni ou bidimensionnel) pour représenter la fissure. La mise à jour de la fissure se fait
simplement en propageant le maillage représentant la fissure.

[Link] Construction de la fissure

La position des éléments unidimensionnels qui représentent la fissure est repérée par
les cinq données suivantes (voir F IG .2.10) :
– Numéro de l’élément coupé iEL ,
– Numéro du côté entrant,
– Position relative par rapport aux nœuds du côté entrant (0 6 αe 6 1),
– Numéro du côté sortant,
– Position relative par rapport aux nœuds du côté sortant (0 6 βs 6 1).

Les côtés de l’éléments sont numérotés de I à IV dans l’ordre de définissions des

noeuds. Les positions relatives des points entrant et sortant de la fissure sont définies
par rapport au premier nœud de chaque côté coupé à l’aide des fractions αe et βs , en
respectant la même convention. Sur la figure 2.10, la position des points entrant et sortant
de la fissure (respectivement A et B) est donnée par :
     
xA x4 x1
 yA  = αe  y4  + (1 − αe )  y1  , (2.42)
zA z4 z1
     
xB x2 x3
 yB  = βs  y2  + (1 − βs )  y3  , (2.43)
zB z2 z3
avec (xi , yi , zi ) les coordonnées du nœud i.

67

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

F IGURE 2.10: Représentation de la fissure sous forme d’un maillage d’éléments

unidimensionnels sur le plan moyen de la coque

La fissure est donc constituée de plusieurs segments dont la position caractérisés par
ces cinq informations, ce qui suffit à complètement la localiser au sein du maillage. Un
exemple simple est donné ci-dessous.

Exemple : Fissure constituée de quatres éléments unidimensionnels placée dans un

maillage de 3x3 quadrangles dont les données sont regroupées dans TAB .2.4 et le résultat
est visible figure 2.11.

TABLEAU 2.4: Données nécessaires pour localiser la fissure.

n˚ de l’élément 1D Élément coupé iEL Côté entrant αe Côté sortant βs

1 4 IV 0.75 II 0.25
2 5 IV 0.5 III 0.5
3 8 I 0.5 II 0.25
4 9 IV 0.75 III 0.75

68

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

F IGURE 2.11: Fissure décrite par TAB .2.4.

Cette description permet également de faciliter le passage de l’élément déformé à

l’élément de référence. Cela est par exemple utile pour étudier la position relative des
nœuds par rapport à la fissure et donc pour leur attribuer leur enrichissement. En effet,
la gestion de la normale et par conséquent la différenciation d’un côté et de l’autre de la
fissure n’est pas chose aisée sur une surface tridimensionnelle courbe. La position relative
des points d’intégration par rapport à la fissure, nécessaire à l’intégration numérique des
forces internes, est également déterminée de la même manière.

Un repère local R =(t1 ,t2 , n) est attribué en chaque pointe de fissure (voir F IG .2.12
avec t1 la direction de la fissure, n la normale à l’élément au point où sort la fissure et t2
un troisième vecteur pour former un repère orthonormé directe.

F IGURE 2.12: Représentation de la fissure sous forme d’un maillage d’éléments

unidimensionnels sur le plan moyen de la coque

69

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

2.2.6 Cas test à fissure fixe

Un cas test à fissure fixe est réalisé dans Europlexus en comparant :
– Une fissure explicitement maillée avec 27x28 éléments de coque classique Q4GS,
– Une fissure dans un maillage composé de 27x27 élément de coque X-FEM XQ4GS.

La structure carrée de côté L = 3m et d’épaisseur e = 0.1m est d’abord soumise à un

chargement en membrane puis à une force ponctuelle hors plan en une de ses extrémité
notée P. Les conditions aux limites imposées sont visibles figure 2.13.

TABLEAU 2.5: Paramètres matériaux utilisés pour la simulation du cas test à

fissure fixe.

Module d’Young E 199 GPa

Coeff. de Poisson ν 0.3
Masse volumique ρ 8000 kg/m3
Limite élastique σy 1990 MPa
Module tangent Et 70 MPa

Ces simulations numériques présentées figure 2.14 affichent des résultats identiques.
Les courbes de déplacement au point P sont présentée figure 2.15 à figure 2.17.

Remarque : Sur la simulation utilisant des éléments de coques XQ4GS, les lèvres de la
fissures semblent liées par les éléments dans lesquels la fissure trouve. Le logiciel de
postraitement Paraview continue de lier les éléments à l’affichage sans savoir qu’ils sont
rompus.

Cet exemple à fissure fixe nous permet de valider les modifications apportées à
l’élément fini de coque pour que celui-ci soit capable de représenter à la fois des
discontinuités en déplacement et celles en rotation. Il s’agit désormais de choisir quand
et comment va avancer la fissure.

70

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

F IGURE 2.13: Schéma du cas test à fissure fixe pour 3x3 éléments coques XQ4GS
(à gauche) et 3x4 éléments coques Q4GS (à droite) - Conditions aux limites.

F IGURE 2.14: Essai de traction hors plan à fissure fixe sur coque X-FEM et coque
classique - Fissure représentée par la méthode X-FEM (à gauche) et de manière
explicite (à droite).

71

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

2.5
Ux FEM
Déplacement (10−3 m)

2 Ux XFEM

1.5

0.5

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Temps (10−2 s)

F IGURE 2.15: Déplacement Ux - Cas test fissure fixe.

Uy FEM
3 Uy XFEM
Déplacement (10−3 m)

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Temps (10−2 s)

F IGURE 2.16: Déplacement Uy - Cas test fissure fixe.

72

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Couplage de la méthode X-FEM au modèle coque de Mindlin-Reissner

20
Uz FEM
Uz XFEM
Déplacement (10−3 m)

15

10

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Temps (10−2 s)

F IGURE 2.17: Déplacement Uz - Cas test fissure fixe.

73

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

2.3 Critère de propagation de fissures en dynamique ex-

plicite
Les aciers utilisés sur les coques de navires ont un comportement qualifié de ductile
à température ambiante et le choix d’un critère de propagation adapté est primordial. Un
des objectifs de cette thèse est d’adapter aux coques un ou plusieurs critères développés
dans la littérature pour permettre de répondre à trois questions :
– Est-ce que la fissure se propage ?
– Dans quelle direction se propage la fissure ?
– A quelle vitesse se propage-t-elle ?

Cette partie énonce tout d’abord deux critères récemment développés pour de la fis-
suration de structures navales. La méthodologie employée pour adapter un de ces critères
à un modèle de coque est ensuite présentée, et enfin la démarche envisagée pour la suite
des travaux est décrite.

2.3.1 Choix d’un critère de propagation

Critère basé sur un modèle élasto-plastique La mécanique de la rupture propose des
critères macroscopiques de propagation de fissure permettant seulement de prédire le
comportement d’une fissure soumis à un mode de fissuration en traction et dans un milieu
fragile. L’idée proposée par Haboussa et al. [HAB 12b] est de développer un critère de
direction de propagation en mode mixte. Il permet ainsi de modéliser la transition d’un
mode de traction à un mode de cisaillement.

Ce modèle utilise une loi de comportement élasto-plastique avec un critère de rupture

basé sur le calcul d’un champ de contrainte équivalent σ̃i j . Celui-ci est évalué à l’aide
d’une pondération sur un demi-disque en avant de la pointe de fissure 2.18 ce qui permet
de s’affranchir de la singularité du champ de contrainte.
Z
r 2
σi j (M)e−( R ) dSM
D
σ̃i j = Z
r 2
, (2.44)
e−( R ) dSM
D

où D est le demi-disque centré sur la pointe de fissure et paramétré par son rayon R
(généralement égal à trois fois la taille d’un élément en pointe de fissure) et r est la
distance entre le front de fissure et le point M.

Ainsi, plus un point de Gauss est proche de la pointe de fissure, plus sa contribution
dans la pondération est importante. De ce champ sont tirées les contraintes principales σ̃I
et σ̃II . Il y a alors propagation si la contrainte principale maximale est supérieure à une
valeur critique σc (ε̇) qui dépend de la vitesse de déformation en pointe :

74

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

F IGURE 2.18: Domaine d’intégration D centré sur la pointe de fissure.


σ̃eq = max(σ̃I , σ̃II ) < σc (ε̇) ȧ = 0 pas de propagation,
(2.45)
σ̃eq = max(σ̃I , σ̃II ) > σc (ε̇) ȧ > 0 propagation.

Critère basé sur un modèle élasto-plastique-endommageable Un second critère dé-

veloppé par J.P. Crété [CRE 13] est basé cette fois-ci sur une loi matériau élasto-plastique
endommageable avec un endommagement de type Gurson modifié [NEE 87], encore
appelé modèle GTN (Gurson-Tvergaard-Needleman). Cela permet entre autre une phase
d’amorçage qui fait la transition entre un endommagement diffus et la formation d’une
fissure. Ce critère est implémenté dans le logiciel industriel Abaqus et utilisé avec un
schéma implicite avec une loi de comportement en déformations planes.

Dans cette approche, une première étape d’identification des plans de localisation
à partir de l’opérateur tangent est réalisée via une analyse de bifurcation. Cela permet
de déterminer deux directions potentielles de propagation 2.19a. Ensuite, en se basant
sur un calcul de champ de déformation équivalent sur une pastille en pointe de fissure,
la direction qui maximise une déformation équivalente est choisie. Pour permettre de
réduire la dépendance à la taille de maille qui apparaît par exemple avec une porosité à
rupture, le critère d’amorçage est basé sur une énergie stockée en pointe 2.19b. Il y a
alors avancée de fissure si le critère est vérifié 2.19c.

Les développements présentés dans cette thèse étendent le premier critère proposé par
Haboussa et al. aux coques. Ce critère a été privilégié car il est plus simple à implémenter.
Par manque de temps l’adaptation au coques du second critère n’a pas été réalisée.

75

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

(a) Détermination de la (b) Vérification du critère sur (c) Calcul de l’avancée de

direction . le demi-disque. fissure.

F IGURE 2.19: Schématisation du modèle de propagation [CRE 13].

2.3.2 Direction de propagation en traction/cisaillement

Quelques notations qui serviront dans la suite du chapitre sont d’abord introduites puis
les direction de propagation établies par Haboussa et al. dans le cas tridimensionnel et ses
simplifications dans le cas bidimensionnel sont présentés. Enfin, l’algorithme de propa-
gation qui permet de basculer automatiquement d’une propagation en mode de traction à
un mode de propagation en cisaillement est synthétisé.

[Link] Notations et changements de bases

Tout d’abord, définissons l’ensemble des repères (représentés F IG .2.20) qui inter-
viennent :
– G =(ex , ey , ez ) définit le repère global du problème,
– R =(t1 ,t2 , n) définit le repère local en pointe de fissure (identique au repère de défi-
nition des modes de fissuration),
– R 0 =(t10 ,t20 , n0 ) définit le repère local en pointe de fissure tourné de l’angle θc =
(t1 ,t10 ) = (t2 ,t20 ), appelé angle de branchement.
– R 00 =(t100 ,t200 , n00 ) est le repère local en pointe de fissure tourné de l’angle θc puis de
l’angle ψc = (t20 ,t200 ) = (n0 , n00 ), ce dernier étant appelé angle de déversement.
Au départ, les contraintes calculées sont exprimées dans le repère local de la pointe
de fissure R . Ses composantes sont alors données par :

 
σ11 σ12 σ13
σi j = σ21 σ22 σ23  . (2.46)
σ31 σ32 σ33

76

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

F IGURE 2.20: Illustration des différents repères utilisés en pointe de fissure.

Les composantes du tenseur des contraintes en pointe de fissure dans le repère R 0 sont
celles du tenseur exprimé dans le repère cylindrique. Ces composantes sont données par :
         
0 K θ 3θ KII θ 3θ
σ11 = σrr = √I 5 cos − cos − √ 5 sin − 3 sin , (2.47)
4 2πr 2 2 4 2πr 2 2
         
0 K θ 3θ KII θ 3θ
σ22 = σθθ = √I 3 cos + cos − √ 5 sin + 3 sin , (2.48)
4 2πr 2 2 4 2πr 2 2
         
0 K θ 3θ KII θ 3θ
σ12 = σrθ = √I sin + sin + √ cos + 3 sin , (2.49)
4 2πr 2 2 4 2πr 2 2
 
0 K θ
σ13 = σrz = √III sin , (2.50)
4 2πr 2
 
0 K θ
σ23 = σθz = √III cos , (2.51)
4 2πr 2
(  0 0

0 ν σ11 + σ22 en déformations planes,
σ33 = σzz = (2.52)
0 en contraintes planes.

Si on calcule maintenant les composantes du tenseur des contraintes dans R 00 , on

obtient en déformation planes les équations 2.53 à 2.58. Ces contraintes sont exprimées
dans le repère en pointe de fissure tourné de θ et de ψ. En contraintes planes, ces équations

77

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

0
sont simplifiées en prenant en compte la relation σ33 = σzz = 0.
00 0
σ11 = σ11 , (2.53)
00 0 2 0 0 2
σ22 = σ33 sin (ψ) + 2σ23 sin (ψ) cos (ψ) + σ22 cos (ψ) , (2.54)
00 0 0 0
σ33 = σ33 cos2 (ψ) − 2σ23 sin (ψ) cos (ψ) + σ22 sin2 (ψ) , (2.55)
00 0 0
σ12 = σ13 sin (ψ) + σ12 cos (ψ) , (2.56)
00
 0 0
 0
σ23 = σ33 − σ22 sin (ψ) cos (ψ) + σ23 2 cos2 (ψ) − 1 ,


(2.57)
00 0 0
σ33 = σ13 cos (ψ) − σ12 sin (ψ) . (2.58)

Enfin, la notion de facteurs d’intensité des contraintes (FIC) normalisés Kni est intro-
duite telle que :
|KI |
KnI = , (2.59)
|KI | + |KII | + |KIII |
|KII |
KnII = , (2.60)
|KI | + |KII | + |KIII |
|KIII |
KnIII = . (2.61)
|KI | + |KII | + |KIII |

Cette notion proposée par Schöllmann et al. [SCH 02] permet d’interpréter les équa-
tions décrivant les directions angulaires dépendant de la mixité tridimensionnelle du
mode de chargement. Afin de visualiser une direction de propagation associée à un tel
chargement mixte, un système de coordonnées barycentrique est utilisé. Haboussa et
al. [HAB 12b] reprennent cette idée en adoptant une approche locale (voir paragraphe
[Link]) compte tenu du caractère non linéaire du matériau en pointe de fissure. Ils re-
lient les FIC aux contraintes équivalentes σ̃i j évaluées sur le demi-disque en pointe de la
manière suivante : 
 KI → σ̃22 ,
KII → σ̃12 , (2.62)
KIII → σ̃23 ,


ce qui conduit à la définition de contraintes équivalentes normalisées σ̃ni :

|σ̃22 |
σ̃n1 = , (2.63)
|σ̃22 | + |σ̃12 | + |σ̃23 |
|σ̃12 |
σ̃n2 = , (2.64)
|σ̃22 | + |σ̃12 | + |σ̃23 |
|σ̃23 |
σ̃n3 = . (2.65)
|σ̃22 | + |σ̃12 | + |σ̃23 |

78

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

Remarque : On rappelle que les tenseurs des contraintes notés σi j sont relatifs à un point
d’intégration, alors que le tenseur des contraintes notés σ̃i j fait référence à une grandeur
équivalente évaluée sur un domaine D. Enfin les contraintes normalisées σ̃ni sont des
grandeurs sacalaires sans dimensions permettant de caractériser la prépondérance d’un
mode par rapport aux autres modes.

[Link] Direction de propagation de fissure tridimensionnelle

Critère de rupture en traction Le critère de direction de propagation proposé par
Erdogan et Sih [ERD 63] (présenté dans le paragraphe [Link]) puis repris par Maigre et
al.[MAI 93] en dynamique est étendu au cas tridimensionnel. Haboussa fait l’hypothèse
que la rupture fragile est régie en statique et en dynamique par le maximum de la
contrainte de traction au voisinage du front de fissure.
00 00
Après avoir montré que σ22 est strictement supérieur à σ33 , l’auteur recherche le
couple d’angles critiques (θc ,Ψc ) qui vérifie :
 00
 ∂σ22 (θ,Ψ)
= 0,
00
∂θ (2.66)
 22∂σ (θ,Ψ)
= 0.∂Ψ

Les développements mènent à une équation polynomiale du troisième ordre en tan θ2
dont la solution n’est pas triviale. Il donne alors les valeurs approchées de θtrac
c et Ψtrac
c
en fonction des FIC normalisés Ki : n

  q 
trac 1
2
θc = 2 sign(KII ) arctan K̂ − K̂ + 8 , (2.67)
4
0
!
2 (θtrac )
1 σ 23 c
Ψtrac
c = sign(KIII ) arctan 0 0 , (2.68)
2 σ22 (θtrac trac
c ) − σ33 (θc )

avec :
1 √

1 + KnI − (1 − KnIII ) p(ν)  ( π − 5ν)
 en déformation plane,
K̂ = où p(ν) = 4 √π
KnII 
 en contrainte plane,
4
(2.69)
0 0
et σ̃i j (θtrac
c ) le tenseur des contraintes équivalent tourné dans le repère R et évalué pour
l’angle θtracc .

Enfin, il substitue les FIC Ki et FIC normalisés Kni par les contraintes équivalentes
σ̃i j et les contraintes équivalentes normalisées σ̃ni à l’aide de la relation 2.62. Il obtient

79

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

finalement :
  q 
1
2
θtrac
c = 2 sign(σ̃12 ) arctan Ŝ − Ŝ + 8 , (2.70)
4
0
!
1 2 σ̃23 (θtrac
c )
Ψtrac
c = sign(σ̃23 ) arctan 0 0 , (2.71)
2 σ̃22 (θtrac trac
c ) − σ̃33 (θc )

avec :
p(ν)
1 + σ̃n1 − 1 − σ̃n3
Ŝ = . (2.72)
σ̃n2

Critère de rupture en cisaillement En présence d’une plasticité non négligeable en

pointe de fissure, Haboussa fait l’hypothèse que la direction de propagation est pilotée
par le maximum de la contrainte équivalente de Von Mises donnée par l’équation 2.73 :
1
q
σVeqM = √ (σ11 − σ22 )2 + (σ22 − σ33 )2 + (σ33 − σ11 )2 + 6(σ212 + σ213 + σ223 ). (2.73)
2
L’auteur recherche le couple d’angles critiques (θcisail
c ,Ψcisail
c ) qui vérifie :
( ∂σV M (θ,Ψ)
eq
= 0,
∂θ
∂σVeqM (θ,Ψ)
(2.74)
∂Ψ = 0.

De la même manière qu’en traction, il exprime les valeurs approchées θcisail

c et Ψcisail
c
n
en fonction des en fonction des FIC normalisés Ki :

– pour KI 6= 0 et KII 6= 0,
KnIII
  
  n
K (1 − )
π 8 3  I 20 
θcisail = sign(KII ) 1+ ν arctan    , (2.75)

c
4 145 5 KnII
0 0
!
1 σ33 (θcisail
c ) − σ22 (θcisail
c )
Ψcisail
c = sign(KIII ) arctan 0 . (2.76)
2 2 σ23 (θcisail
c )

– pour KI = 0 et KII = 0,

θcisail
c = sign(KII ) 0.536, (2.77)
0 0
!
1 σ33 (θcisail
c ) − σ22 (θcisail
c )
Ψcisail
c = sign(KIII ) arctan 0 cisail
. (2.78)
2 2 σ23 (θc )

Les FIC Ki et FIC normalisés Kni sont à nouveau substitués par les contraintes
équivalentes σ̃i j et les contraintes équivalentes normalisées σ̃ni pour donner :

80

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

– pour KI 6= 0 et KII 6= 0,
σ̃n3
  
  n
σ̃ (1 − )
π 8 3  1 20 
θcisail = sign(σ̃12 ) 1+ ν arctan    , (2.79)

c n
4 145 5 σ̃2
0 0
!
1 σ̃33 (θcisail
c ) − σ̃22 (θcisail
c )
Ψcisail
c = sign(σ̃23 ) arctan 0 cisail
. (2.80)
2 2 σ̃23 (θc )

– pour KI = 0 et KII = 0,

θcisail
c = sign(σ̃12 ) 0.536, (2.81)
0 0
!
1 σ̃33 (θcisail
c ) − σ̃22 (θcisail
c )
Ψcisail
c = sign(σ̃23 ) arctan 0 cisail
. (2.82)
2 2 σ̃23 (θc )

Transition entre les deux critères Les différentes directions de propagation possibles
en traction (θtrac ) ou en cisaillement (θcisaillement ) sont calculées à partir du champ de
contrainte équivalent σ̃i j . La détermination de la direction de propagation est réalisée en
estimant l’intensité de la déformation plastique équivalente en pointe de fissure. Pour ce
faire, un tenseur des déformations équivalent en pointe de fissure ε̃i j est calculé de la
même manière que pour la contrainte équivalente :
Z
r 2
εi j (M)e−( R ) dSM
D
ε̃i j = Z . (2.83)
−( Rr )2
e dSM
D

Une déformation équivalente au sens de Mazars [MAZ 84] notée ε̃eq est calculée.
Celle-ci s’appuie s’appuie sur les deformations principales de traction et est définie de la
manière suivante : p
ε̃eq = < ε̃I >2 + < ε̃II >2 . (2.84)
Si cette déformation équivalente est inférieure à un certain seuil ε̃trac , la direction
sera pilotée par le mode de rupture en traction (voir F IG .2.21). Si au contraire cette
déformation est supérieure à ε̃cisail , la direction sera donnée par le cisaillement. Entre ces
deux seuils, une loi des mélanges entre les deux directions est appliquée pour calculer un
angle intermédiaire θtrans donné par l’équation 2.85.

θtrans = Mel θcisail

c + (1 − Mel ) θtrac
c , (2.85)
où :
ε̃eq − ε̃trac
Mel = .
ε̃cisail − ε̃trac
L’identification des paramètres ε̃trac , ε̃cisail et σc (ε̇) est basée sur une méthode dite
inverse. Il faut pour cela disposer d’une expérience pour laquelle l’historique du front

81

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

F IGURE 2.21: Influence de la déformation équivalente pour choisir la direction de

propagation.

de fissure, autrement dit la position de la pointe en chaque instant est connue. Cela
s’obtient en général à partir de l’exploitation des images en caméra rapide. Ensuite, cette
expérience est simulée en imposant l’ouverture spatio-temporelle expérimentale. Avec
les résultats des simulations, on peut alors prétendre connaître les champs mécaniques, à
chaque instant t, qui ont permis de faire propager ou non la fissure expérimentale. Enfin,
on peut identifier les différents paramètres à partir de l’exploitation de ces champs. En
pratique, cette identification n’est pas simple car elle nécessite un recalage très précis
de l’essai et de la simulation sans quoi les résultats sont rapidement erronés. A noter
également que la connaissance de la loi de comportement viscoplastique du matériau
pour la température de l’essai est nécessaire.

Une fois que les paramètres ont été identifiés, le critère peut être utilisé dans les si-
mulations pour réaliser des calculs prédictifs. Cependant il faut veiller à toujours garder
une taille de maille comparable à celle utilisée lors de l’identification. Les paramètres dé-
pendent en effet de la taille de la zone de régularisation D qui est calculée sur un nombre
prédéfini d’éléments. Haboussa a montré que pour pouvoir utiliser les paramètres identi-
fiés sur un premier maillage dans une autre simulation, il faut veiller à conserver la même
taille physique de zone D.

[Link] Simplifications dans le cas bidimensionnel

Dans le cas bidimensionnel l’angle de déversement Ψ n’intervient plus dans le pro-

blème et les équations sont simplifiées. Pour la rupture en mode de traction, Haboussa
reprend la solution analytique proposée par Erdogan et Sih [ERD 63]. Pour la rupture en
mode de cisaillement, il cherche la direction qui maximise la contrainte de cisaillement et

82

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

obtient cette fois-ci une solution analytique donnée par l’équation 2.86.
 h i
x2
 r   
2x 1 +
θtrac  3 R (x) cos  1 arccos 
anal = 2 arctan 2
9
 − 2π  + x  , (2.86)
6 3 R (x) 3 3

avec : r
σ̃22 343 49 2 14 4 4 6
x= et R (x) = + x + x + x .
σ̃12 6 3 9 81
Une solution semi-analytique approchée est proposée afin de simplifier l’implémenta-
tion de l’angle en cisaillement. Finalement, il obtient :
  s 
 2
1 σ̃22 σ̃22 
θtrac = 2 arctan   − sign(σ̃12 ) 8 + , (2.87)
4 σ̃12 σ̃12
4 √ σ̃22
 
cisail π
θ = sign(σ̃12 ) arctan π . (2.88)
4 19 σ̃12

Le même principe de transition que celui présenté dans le cas tridimensionnel est
utilisé.

[Link] Algorithme de transition traction/cisaillement

La transition d’un mode de rupture à un autre se fait en étudiant l’intensité de la défor-
mation plastique en pointe de fissure. Haboussa l’évalue à l’aide d’une déformation équi-
valente présentée dans la partie [Link]. Une synthèse du critère en traction/cisaillement
est proposé dans l’algorithme 1.

2.3.3 Incrément de propagation dans un pas de temps

La vitesse de propagation ȧ est calculée à partir de la relation de Kanninen et
Popelar [KAN 85] présentée équation 1.29, en remplaçant le facteur d’intensité des
contraintes circonférentielles Kθθ par la contrainte équivalente σ̃eq et la ténacité KIc par
une contrainte critique σc (ε̇) dépendant de la vitesse de déformation :
 
σc (ε̇)
ȧ = cR 1 − . (2.89)
σ̃eq
Ainsi, connaissant le pas de temps, on remonte à l’avancée de fissure ∆a.

Toutes les équations décrivant ce modèle de propagation de fissure sont basées sur le
concept de l’approche locale en pointe de fissure et dérivent du modèle énergétique. Ainsi
la contrainte critique σIc peut être considérée symboliquement comme une mesure locale
du taux de restitution de l’énergie critique GIc .

83

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

Algorithme 1: Algorithme de transition traction-cisaillement 2D

Entrées : σ̃i j , ε̃i j ,σIc , cr , ε̃cisail , εtract .
Sorties : ȧ, θc .
Calcul des tenseurs équivalents σ̃i j et ε̃i j sur le domaine D :
R R
D σRi j (M)w(r)dD D εRi j (M)w(r)dD
σ̃i j = , ε̃i j =
D w(r)dD D w(r)dD

Évaluation des critères :

p
σ̃eq = max(σ̃I , σ̃II , 0), ε̃eq = {ε̃I }+ , {ε̃II }+ .

si (σ̃eq ≥ σIc ) alors

Déterminer
 la vitesse de propagation de la fissure :
σIc
ȧ = cr 1 − ; . voir Eq. 2.89
σ̃eq
Calculer les directions
  de propagations
: 
 2
1 σ̃22 σ̃22
θtrac = 2 arctan   − sign(σ̃12 ) 8 +  ; . voir Eq. 2.87
4 σ̃12 σ̃12
4 √ σ̃22
 
π
θcisail = sign(σ̃12 ) arctan π ; . voir Eq. 2.88
4 19 σ̃12
θtrans = Mel θtrac cisail ;
c + (1 − Mel ) θc . voir Eq. 2.85

Déterminer la direction de propagation :


trac si ε̃ ≤ ε̃
θ
 eq trac ,
θc = θ cisail si ε̃eq ≥ ε̃cisail ,

 trans
θ sinon.
sinon
ȧ = 0, θc = 0.
fin

84

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

2.3.4 Extension du critère pour les coques

Maintenant que le critère d’initiation et de direction proposé dans les travaux
d’Haboussa [HAB 12b] a été présenté, il s’agit de décrire les différentes modifications
nécessaires à son utilisation avec des éléments finis de coque.

La première différence vient de la formulation des éléments finis de coques : le tenseur

des contraintes est exprimé dans le repère local noté L propre à chacun des points d’in-
tégration (voir F IG .2.22). L’évaluation d’un champ de contraintes équivalent en pointe
de fissure nécessite donc une opération de changement de base pour chaque point d’in-
tégration présent dans D afin d’ajouter la contribution du tenseur associé pondéré par la
fonction poids ω.

F IGURE 2.22: Mise en évidence du repère local en pointe de fissure R et des

repères locaux Li associés aux différents points d’intégration.

La deuxième étape est de redéfinir la forme du domaine D utilisé pour calculer le

tenseur des contraintes moyenné en pointe de fissure. Dans le cas bidimensionnel, Ha-
boussa [HAB 12b] utilise un demi-disque centré sur la pointe de fissure située devant
celle-ci. Pour l’extension au cas tridimensionnel, Pelée de Saint Maurice [SAI 14] étudie
l’influence de différentes formes de domaines (demi-cube, demi-cylindre et demi-sphère)
sur le résultat. L’auteur a entre autres montré qu’il est plus judicieux d’utiliser une demi-
sphère dans ce cas là.

Dans le cas des coques, bien que la géométrie de l’élément soit représentée par une
surface, le champ de contrainte peut être évalué et varier dans l’épaisseur de celle-ci.
L’élément de coque enrichi X-Q4GS construit précédemment possède cinq points
d’intégration à travers l’épaisseur qui ne sont activé qu’en cas de non linéarité matériau.
Ainsi, pour rajouter de l’information par rapport au problème plan, l’idée n’est plus
de sommer les contributions sur un demi-disque en avant uniquement, mais sur cinq
demi-disques comme on peut le voir figure 2.23.

85

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

F IGURE 2.23: Calcul de la contrainte équivalente pour les coques - Forme du

domaine D pour coque plane.

Lorsque la coque est courbée, le domaine D ne correspond en réalité pas exactement

à un demi disque mais à une géométrie s’en rapprochant. L’algorithme proposé cherche
dans un premier temps les points de Gauss situés sur le plan moyen de coque et inclus
dans une demi-sphère devant la pointe, puis remonte aux contraintes disponibles aux
points de Simpson dans l’épaisseur pour les points des Gauss concernés. La figure 2.23
présente donc la forme du domaine D dans le cas d’une coque plane pour des matériaux
linéaires et non linéaires.

Plusieurs pistes peuvent alors être envisagées pour adapter le critère à l’utilisation des
différentes couches :
– Calculer des champs équivalents sur chacun des cinq disques séparément, évaluer
chacune des contraintes équivalentes associées et propager la fissure dès lors que le
critère est atteint sur une des couches, sur au moins trois couches, ou bien sur toutes
les couches.
Z
σi j (M)w(r)dD
Dk
σ̃i j (couchek ) = Z , Pour k = {1..5}, (2.90)
w(r)dD
Dk

– Calculer un champ équivalent intégré dans l’épaisseur en utilisant la pondération

classique des points de Simpson, évaluer la contrainte équivalente issue de ce

86

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

champ et utiliser le critère classique.

Z Z h/2
σi j (M)w(r)dzdD
D −h/2
σ̃i j = Z Z h/2 (2.91)
w(r)dzdD
D −h/2

Dans nos travaux, c’est cette deuxième méthode qui a été implémentée.

2.3.5 Moyennes temporelles sur les avancées du front de fissure

A chaque pas de temps ∆t, le champ de contrainte équivalent σ̃i j est évalué sur le
domaine D intégré dans l’épaisseur. Si le critère est vérifié, une avancée ∆a dans une
direction θ est calculée avec les équations présentées dans 57. Cependant, quelques pré-
cautions doivent être prises avant de mettre à jour la fissure.
Dans la partie [Link], il a été choisi de ne garder que les enrichissements de type saut
et de ne pas utiliser les enrichissements asymptotiques en pointe. Par conséquent, la pointe
de fissure ne peut pas être localisée à l’intérieur d’un élément. En d’autres termes, si le
critère produit une avancée ∆a1 plus petite que l’élément pour un pas de temps ∆t1 , tout
se passe mathématiquement comme si l’élément était entièrement coupé (voir F IG .2.24).
Cela peut à la fois conduire à une dissipation trop importante dans la fissuration, mais
également à une vitesse de propagation de fissure supérieure à la célérité des ondes.

F IGURE 2.24: Mise à jour de la fissure pour une avancée suffisante.

Le pas de temps ne pouvant pas être augmenté, il est nécessaire d’attendre que la
longueur d’avancée produite par le critère soit de l’ordre de grandeur de la taille d’un
élément notée Lcar .

87

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

Vient alors la question du choix de la direction dans laquelle propager la fissure. En

effet, d’un pas de temps à l’autre, le champ de contrainte peut fluctuer entrainant égale-
ment une variation de l’angle de propagation. Prendre l’angle obtenu au dernier pas de
temps conduirait à ne pas tenir compte de l’histoire du chargement. La première idée est
de réaliser une moyenne sur les angles de propagation obtenus pendant la phase où le cri-
tère est atteint, mais la moyenne temporelle des champs équivalents qu’on utilise ensuite
pour calculer l’angle de propagation semble donner de meilleurs résultats.
Ainsi, on choisit d’effecter cette moyenne temporelle sur les champs équivalents qu’on
moyT
note σ̃i j et donnée par :
 
moyT ∆t moyT ∆t
σ̃i j = 1− σ̃i j + σ̃i j , (2.92)
Tcumul Tcumul

où Tcumul est un compteur qui cumule le temps à partir du moment où le critère est atteint et
tant que celui-ci est vérifié, jusqu’à ce qu’une avancée suffisante soit atteinte (F IG .2.25).
Cette avancée est donnée par l’équation 2.93.

∆a = ȧ Tcumul , (2.93)
où ȧ est la vitesse calculée avec la relation de Kanninen (équation 2.89) en utilisant la
contrainte équivalente issue de la moyenne temporelle des champs équivalents.

F IGURE 2.25: Mise à jour de la fissure pour une avancée suffisante.

Le compteur Tcumul est réinitialisé à chaque propagation effectuée mais également dès
que le critère n’est plus vérifié. Cela permet par exemple d’éviter la propagation artificielle
d’une fissure dans le cas où le critère de propagation n’est dépassé que sur quelques pas de
temps. Les modifications liées aux moyennes temporelles sur le critères sont présentées
dans l’algorithme 2.

88

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

Algorithme 2: Modifications liées au moyennes temporelles sur le critère.

si (σ̃eq < σIc ) alors
Pas de propagation :
Propagation = False ;
ȧ = 0 ;
∆a = 0 ;
moyT
σ̃i j = 0 ;
Tcumul = 0 ;
sinon
Incrémentation du compteur et moyennes temporelles :
Tcumul = Tcumul + ∆t ;
moyT
Moyenne temporelle du tenseur équivalent σ̃i j ; . voir Eq. 2.92
moyT moyT moyT
Evaluation du critère σ̃eq = max(σ̃I , σ̃II , 0) ; . voir Eq. 2.45
σIc
Vitesse de propagation : ȧ = cr (1 − moyT ) ; . voir Eq. 2.89
σ̃eq
Calcul de l’avancée : ∆a = ȧ Tcumul ; . voir Eq. 2.93
si ∆a > Lcar alors
Propagation = True ;
fin
fin

2.3.6 Mise en œuvre numérique de la propagation

A partir des différents champs mécaniques, le critère permet de donner une direction
de propagation θ et une avancée de fissure ∆a cohérente avec la taille des éléments définie
précédemment. Il faut alors mettre à jour la géométrie de la fissure puis activer et évaluer
les enrichissements (voir Schéma 2.26). La méthodologie employée pour ces deux étapes
est présentée dans les paragraphes [Link] et [Link].

F IGURE 2.26: Schématisation des différentes étapes nécessaire au calcul d’une

avancée de fissure.

[Link] Mise à jour de la géométrie de la fissure

Comme décrit dans la partie [Link], la fissure est discrétisée par un maillage
d’éléments unidimensionnels évoluant sur la surface moyenne de la coque, de bord
d’élément à bord d’élément. Bien que plus simple que l’actualisation de level-set, cette

89

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

étape n’est pas pour autant triviale. En effet, les éléments de coques peuvent être gauches
et le trajet de fissure doit alors l’être lui aussi. Pour contourner cette difficulté, la position
de la nouvelle pointe de fissure se fait en passant par l’élément de référence qui est plan.
Les différentes étapes de ce calcul sont présentées ci-dessous.

Pour illustrer ce cas de figure, prenons une fissure dont la pointe est notée A, munie
d’un repère local orthonormé direct R = (A ;~t1 ,~t2 ,~n) et devant propager dans un élément
courbé qui n’est pas dans le même plan que le précédent (voir F IG .2.27). Le critère donne
un angle de propagation θ dans le plan (A ;~t1 ,~t2 ), repéré par rapport à l’axe (A,~t1 ) et avec
une avancée de ∆a.

F IGURE 2.27: Schéma d’une fissure devant se propager dans un élément courbe.

Comme on peut le voir figure 2.28, le point M engendré par cette avancée n’est pas
sur la surface de la coque. Il est alors nécessaire de trouver une manière d’obtenir l’inter-
section avec un des bords de l’élément suivant, synonyme de nouvelle pointe de fissure.
Cet élément étant courbe, deux obstacles se dressent :
– La normale à l’élément est différente en tout point or il faut en choisir une pour
réaliser la projection,
– La surface de projection n’est pas plane.

L’idée retenue est d’utiliser le plan qui contient le côté d’élément sur lequel se trouve
la pointe de fissure et le premier côté adjacent dans le sens trigonométrique. Sur la figure
2.29, ce plan est appelé noté P = (N1 ; e~1 , e~2 ). Quant à la normale de projection, c’est le
vecteur ~n du repère local en pointe qui est choisi.

L’étape suivante est de projeter un point P appartenant au segment de fissure [AM] sur
le plan P décrit précédemment suivant la normale~n. On obtient alors le point P0 qui a pour
coordonnées (α, β, 0) dans le repère orthonormé direct (N1 ; e~1 , e~2 , e~3 ). La droite (AP0 ) ne
peut pas directement être prise comme direction de la fissure car elle ne s’intersecte pas

90

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

F IGURE 2.28: Visualisation de l’angle et de l’avancée de fissure donnés par le

critère.

F IGURE 2.29: Plan choisi pour la projection dans l’étape de propagation noté P .

91

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

avec un des côtés de l’élément (côté N2 N3 sur la figure 2.30).

F IGURE 2.30: Projection d’un point du trajet de la fissure sur un plan passant par
trois points de l’élément.

Pour s’affranchir de ce problème, une solution consiste à utiliser la bijection qui existe
entre l’élément déformé et l’élément de référence. Ce dernier étant plan, le calcul d’in-
tersection y est trivial. L’image du point P0 dans l’élément de référence notée Pre 0
f est
obtenue via les équations 2.94 et 2.95 qui traduisent simplement la conservation de la
position relative du point par rapport à un noeud dans l’élément (voir F IG .2.31). Il en est
de même pour le point Are f image de la pointe de fissure A.
α αre f
−−−→ = −→ , (2.94)
kN1 N2 k kI1 I2 k
β βre f
−−−→ = −→ . (2.95)
kN1 N4 k kI1 I4 k
0 ) et un des côtés de cet élément donne le
L’intersection entre la demi-droite [Are f Pre f
point Bre f . La positon du point sur le maillage déformé dont l’image est Bre f dans l’élé-
ment de référence est obtenu à l’aide des fonctions d’interpolation géométrique classique :

hxB i = ∑ Ni (ξBre f , ηBre f )hXi i, (2.96)

i=1,4

hxB i = hxB yB zB i ; hXi i = hXi Yi Zi i, (2.97)

où xB , yB , zB sont les coordonnées cartésiennes de la nouvelle pointe de fissure notée B,
Xi , Yi , Zi sont les coordonnées cartésiennes du nœud i, et Ni (ξ, η) les fonctions de forme

92

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

(a) Position de l’image de P0 . (b) [Are f Pre0 f ] prolongé pour obtenir Bre f .

F IGURE 2.31: Passage dans l’élément de référence pour l’étape de propagation.

F IGURE 2.32: Calcul de la position de la pointe de fissure sur le maillage déformé.

93

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

de l’élément bilinéaire.

La méthodologie qui vient d’être présentée permet donc de mettre à jour la géométrie
de la fissure à partir d’une avancée ∆a et d’un angle θ calculés par le critère. Il convient
de remarquer que l’avancée produite notée ∆ae f f n’est pas rigoureusement la même que
cette calculée par le critère du fait de la géométrie de la coque. On observe en effet sur
l’exemple que le point M calculé (voir F IG .2.28) est différent du point B produit (voir
F IG .2.32). Cependant, la courbure d’un élément tout comme le changement d’orientation
d’un élément à l’autre sont généralement faibles, ce qui rend les approximations accep-
tables.

Il peut également arriver que l’avancée ∆ae f f produite dans l’élément soit plus faible
que celle prévue par le critère (qui est de l’ordre de la taille d’un élément). Cela peut par
exemple se produire lorsque la pointe se trouve initialement proche d’une extrémité du
côté entrant et que la fissure vient intersecter le côté adjacent (voir F IG .2.33). Dans ce
cas, en fonction de l’avancée qui reste à produire pour être cohérent avec le critère, la
fissure propage à nouveau ou non dans l’élément suivant dans la même direction que celle
initialement prévue par le critère. Dans nos travaux, si après une propagation la longueur
qui reste à propager est supérieure à 0.3 Lcar , alors on choisit de couper l’élément suivant.
Cette procédure est décrite dans l’algorithme 3.
Enfin, dans le cas où la fissure propage à nouveau dans le même élément, on décide
de décrire la fissure en prenant le premier point d’entrée et le point de sortie de l’élément
comme on peut le voir figure 2.34.

F IGURE 2.33: Avancée supplémentaire lorsque l’avancée effective ∆ae f f est trop
faible par rapport à ce que prévoir le critère.

94

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Critère de propagation de fissures en dynamique explicite

F IGURE 2.34: Modification de la fissure lorsqu’elle propage dans le même

élément.

Algorithme 3: Propagation pour produire une avancée cohérente avec le critère

Entrées : ∆a, Propagation.
Résultat : Mise à jour de la géométrie de la fissure.
si (Propagation == True) alors
moyT
σ̃i j = 0 ;
Tcumul = 0 ;
tant que (∆a > 0.3 Lcar ) faire
Propager dans l’élément suivant ;
Avancée effective de ∆ae f f ;
∆a = ∆a − ∆ae f f ;
fin
Propagation = False ;
sinon
fin

95

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

[Link] Mise à jour de l’enrichissement

Dès lors que la géométrie de la fissure est établie, la dernière étape à réaliser est de
mettre à jour l’enrichissement. Cette question a déjà été en partie évoquée dans le pa-
ragraphe [Link] concernant la construction de la fissure. La position relative d’un point
d’intégration ou d’un noeud d’un élément par rapport à la fissure est obtenue dans l’élé-
ment de référence. Si on note I un point quelconque de l’élément coupé par la fissure AB,
Ire f le point et Are f Bre f le segment correspondant dans l’élément de référence et →−v un
−−−−−→
vecteur orthogonal à Are f Bre f , l’enrichissement est donné par :
−−−−→ −
+1 si Are f Ire f .→
(
v > 0,
H(xI , yI ) = H(ξIre f , ηIre f ) = −−−−→ → − (2.98)
−1 si Are f Ire f . v < 0.

F IGURE 2.35: Détermination de l’enrichissement à partir de l’élément de référence.

96

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Simulation d’un exemple type issu de la littérature

2.4 Simulation d’un exemple type issu de la littérature

Il s’agit désormais de simuler un essai de propagation proposé dans la littérature afin
de tester l’efficacité du critère après les modifications apportées pour une utilisation dans
les coques. On propose ici de réaliser le cas test de Kalthoff et Winkler [KAL 88] qui ne
sollicite la structure qu’en mode I et mode II. Des cas tests faisant intervenir un charge-
ment en mode III sont présentés dans le C HAP.4.

Description L’expérience réalisée par Kalthoff et Winkler [KAL 88] consiste à impac-
ter une éprouvette rectangulaire possédant deux entailles initiales identiques et parallèles
au moyen d’un projectile cylindrique (F IG .2.36). Suivant la vitesse du chargement, la
nuance du matériau ou la température, différents modes de ruptures apparaissent. Pour
des faibles vitesses, l’expérience montre une rupture dirigée par la traction (mode I) ca-
ractérisée par un angle de propagation moyen d’environ 70◦ par rapport à la direction
initiale de l’entaille. Pour des vitesses d’impact élevées, la fissure se propage en mode
II suivant une bande de cisaillement adiabatique et s’initie avec un angle quasi nul et
propage avec un angle de −10◦ . Ces résultats sont visibles figure 2.37.

F IGURE 2.36: Schéma de principe de l’expérience de Kalthoff et Winkler et

observation de l’intensification du champ de contrainte pour une éprouvette
métallique ou transparente - Schéma tiré de [KAL 00].

Données Seule la moitié de l’éprouvette est modélisée étant donné la symétrie du

problème. L’impact est modélisé par une vitesse V0 imposée (et maintenue) sur le côté

97

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

DAMAGE PATHS

DAMAGE SURFACES

LOW RATES HIGH RATES

F IGURE 2.37: Fissurations par traction et par bande de cisaillement observées sur
de l’acier X2 NiCoMo 1895 (encore appelé acier Maraging 18Ni1900) - Images
tirées de [KAL 00].

gauche de l’éprouvette qui contient une préfissure de longueur L = 50mm. L’épaisseur e

de l’éprouvette est de 16mm. Deux cas de calcul sont présentés : un premier cas d’impact
avec une vitesse de V0 = 16m/s et un second cas pour un impact de V0 = 32m/s. Les
conditions aux limites sont visibles figure 2.38.
Le maillage est composé de 100 × 100 éléments coque X-Q4GS dont la taille est de 1mm
et est présenté figure 2.39.

Le matériau utilisé est un acier Maraging 18Ni1900 avec une loi élastoplastique à
écrouissage isotrope dont les paramètres sont donnés dans le tableau 2.6. Concernant le
critère, les mêmes paramètres que ceux utilisés dans les travaux de Haboussa et al. sont
choisis (voir TAB .2.7).

98

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Simulation d’un exemple type issu de la littérature

F IGURE 2.38: Géométrie et conditions aux limites de la simulation de Kalthoff et

Winkler.

F IGURE 2.39: Maillage pour la simulation de l’essai de Kalthoff et Winkler.

TABLEAU 2.6: Paramètres matériaux de l’acier Maraging utilisés pour

l’expérience de Kalthoff.

Module d’Young E 200 GPa

Coeff. de Poisson ν 0.3
Masse volumique ρ 7830 kg/m3
Célérité des ondes de Rayleigh cR 2800 m/s
Limite élastique σy 2000 MPa
Module tangent Et 1600 MPa

99

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

TABLEAU 2.7: Paramètres du modèle de transition traction/cisaillement de l’acier

Maraging.

Contrainte critique σIc 250 MPa

Seuil en traction ˜
εhoop 8.00 ×10−4
Seuil en cisaillement ˜
εshear 8.85 × 10−4
Longueur caractéristique des éléments Lcar 0.001 m
Rayon du demi disque D R 0.003 m

Résultats Les trajets de fissure obtenus sont affichés figure 2.40 et l’évolution de la
contrainte de von Mises au cours des simulations est donnée figures 2.41 et 2.42.
Pour V0 = 16m/s, on retrouve le trajet caractéristique de la rupture fragile en mode I
obtenu expérimentalement sur cette éprouvette, avec un angle de propagation proche de
70◦ . Ensuite, en modifiant uniquement la vitesse d’impact à V0 = 32m/s, on obtient cette
fois-ci une propagation pilotée par le cisaillement pour une fissure sollicitée en mode II,
avec là encore un trajet de fissure proche de celui obtenu expérimentalement.

On peut remarquer que les résultats ne sont pas rigoureusement identiques à ceux
présentés dans les travaux de Haboussa [HAB 12a]. Cela vient du fait que la loi de com-
portement est ici formulée en contraintes planes ce qui a notamment pour conséquence
de modifier la célérité des ondes mécaniques.

Le cas en cisaillement permet une validation du critère de Haboussa et al.. A noter

que nous avonsconstaté lors de la phase initiale de propagation (i.e. lorsque l’angle
de propagation α ' 0◦ ) qu’un écart même faible entre la solution analytique et semi-
analytique avait un impact important sur le résultat obtenu. Les calculs présentés ici ont
donc été effectués avec l’angle analytique (voir équation 2.86). Nous préconisions son
usage plutôt que celui du critère simplifié.

100

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Simulation d’un exemple type issu de la littérature

10
V0 = 16 m/s
9 V0 = 32 m/s
8

6
Y (cm)

0
0 1 2 3 4 5 6 7 8 9 10
X (cm)

F IGURE 2.40: Simulation de Kalthoff et Winkler - Trajets de fissure pour

V0 = 16m/s et V0 = 32m/s.

101

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

F IGURE 2.41: V0 = 16m/s - Évolution de la contrainte de von Mises sur la

simulation de l’expérience de Kalthoff et Winkler pour t=16 µs, t=24 µs, t=28 µs,
t=36 µs, t=44 µs et t=52 µs.

102

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 Simulation d’un exemple type issu de la littérature

F IGURE 2.42: V0 = 32m/s - Évolution de la contrainte de von Mises sur la

simulation de l’expérience de Kalthoff et Winkler pour t=16 µs, t=24 µs, t=28 µs,
t=32 µs, t=36 µs et t=40 µs.

103

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés
 2. Extension de l’approche X-FEM en dynamique rapide pour la propagation de fissure
dans des coques de type Mindlin-Reissner

2.5 Conclusion
Dans ce chapitre nous avons présenté les développements et les choix réalisés pour
simuler la propagation d’une fissure dans une structure fine.

Après avoir rappelé quelques propriétés de l’élément fini de coque épaisse choisi, sa
formulation a été étendue pour que cet élément soit utilisable avec méthode X-FEM. Les
masses associées aux degrés de liberté enrichis sont prises égales à celles associées aux
degrés de liberté classiques et le pas de temps critique est deux fois plus petit que celui
obtenu avec des éléments finis classiques.

Nous avons ensuite fait le choix de discrétiser la fissure avec un maillage unidi-
mensionnel évoluant sur la surface moyenne de la coque. Cela permet de se passer de
l’utilisation des level-sets et simplifier l’étape de propagation. Une stratégie a été mise en
place afin de mettre à jour la géométrie de la fissure avec la description choisie.

Par la suite, un critère de propagation basé sur l’état de contrainte en pointe de fissure
proposé par Haboussa et al. a été repris puis adapté à l’élément fini de coque développé.
La prise en compte du mode III de sollicitation a notamment été ajouté dans le calcul de
la direction de propagation. Ce critère simplifié permet à la fois de simuler la rupture
en traction ou en cisaillement avec une transition automatique entre les deux modes de
propagation. Le calcul de la contrainte équivalente pour déterminer s’il y a propagation
ou non se fait sur cinq demi-disques superposés dans l’épaisseur, dans le voisinage de la
fissure.

Enfin, la simulation d’un exemple type issu de la littérature est présentée et les résul-
tats affichés correspondent aux résultats expérimentaux.

104

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

© [Y. Jan], [2016], INSA Lyon, tous droits réservés