Modélisation hyperélastique des composites
Modélisation hyperélastique des composites
Thèse
Pour obtenir
Le grade de docteur
Formation doctorale
Ecole doctorale Mécanique, Energétique, Génie Civil, Acoustique (MEGA)
Spécialité : Mécanique, génie mécanique, génie civil
Par
Adrien Charmetant
Jury
P. Boisse Professeur (INSA Lyon) Directeur
M. Brieu Professeur (Ecole Centrale Lille) Rapporteur
P. Ladeveze Professeur (ENS Cachan) Examinateur
S. Lomov Professeur (Katholiek Universiteit Leuven) Rapporteur
A. Madeo Maître de conférences (INSA Lyon) Examinateur
E. Maire Directeur de recherches CNRS (INSA Lyon) Invité
S. Otin Ingénieur de recherches (Snecma) Examinateur
F. Sidoroff Professeur émérite (Ecole Centrale Lyon) Examinateur
E. Vidal-Sallé Maître de conférences (INSA Lyon) Examinateur
Résumé
La simulation des procédés de mise en forme des composites à renforts tissés de type RTM est un enjeu
majeur pour les industries de pointe mettant en œuvre ce type de matériaux. Au cours de ces procédés, la
préforme tissée est souvent soumise à des déformations importantes. La connaissance et la simulation du
comportement mécanique de la préforme à l’échelle macroscopique et à l’échelle mésoscopique s’avère
souvent nécessaire pour optimiser la phase de conception de pièces composites formées par de tels
procédés.
Une analyse du comportement mésoscopique des préformes tissées de composites est d’abord
proposée. Une loi de comportement hyperélastique isotrope transverse est développée, permettant de
décrire le comportement mécanique de chacun des modes de déformation de la mèche : élongation dans
la direction des fibres, compaction et distorsion dans le plan d’isotropie de la mèche, cisaillement le long
des fibres. Une méthodologie est proposée pour identifier les paramètres de cette loi de comportement à
l’aide d’essais sur la mèche et sur le tissu, et une validation par comparaison avec des essais
expérimentaux est présentée.
Une analyse du comportement macroscopique des renforts interlocks est ensuite proposée : une loi
de comportement hyperélastique orthotrope est développée et implémentée. Cette loi, extension de la loi
de comportement pour la mèche, est également basée sur une description phénoménologique des modes
de déformation de la préforme. Une méthode d’identification des paramètres de cette loi de
comportement est mise en œuvre, utilisant des essais expérimentaux classiques dans le contexte des
renforts tissés (tension uniaxiale, compression, bias extension test, flexion). Cette seconde loi de
comportement est validée par comparaison avec des essais de flexion et d’emboutissage hémisphérique.
Mots-Clés
Renforts tissés, composites, mise en forme, propriétés mécaniques, analyse mésoscopique, analyse
macroscopique, lois de comportements hyperélastiques, grandes transformations, méthode par éléments
finis
Keywords: Textile composites, woven reinforcements, resin transfer moulding (RTM), mechanical
properties, mesoscopic analysis, macroscopic analysis, hyperelastic behaviour laws, finite
transformations, finite element analysis (FEA)
Ce premier chapitre est une introduction à l’étude des renforts tissés de composites. Après une
présentation générale des matériaux composites, les renforts textiles de composites seront détaillés.
Dans ce manuscrit, nous serons amenés à nous intéresser au comportement mécanique de ces
matériaux à l’échelle mésoscopique et à l’échelle macroscopique. Les différentes étapes de fabrication
d’un renfort tissé ayant une influence directe sur sa structure et donc sur son comportement
mécanique à chacune de ces échelles, le procédé de fabrication du renfort tissé sera décrit
succinctement, de la fabrication de la fibre au tissage de la pièce finale.
Le comportement mécanique des renforts tissés aux différentes échelles sera ensuite décrit, et les
essais permettant de caractériser les différents modes de déformation du renfort tissé seront présentés.
A l’issue de cette partie, une bonne compréhension des phénomènes et des mécanismes pouvant
apparaître au sein du renfort sera acquise, qui permettra de donner un cadre aux modèles et aux lois
de comportement qui seront développés par la suite.
Finalement, les modèles et lois de comportement permettant à l’heure actuelle de simuler le
comportement mécanique des renforts tissés de composites aux différentes échelles seront présentés,
avec leurs avantages et leurs limites.
Sommaire
Renforts de composites
Les renforts sont généralement classés en fonction de leur géométrie. On distingue notamment :
• les renforts particulaires : il s’agit généralement d’inclusions granulaires (craie), lamellaire
(talc, mica) ou aciculaire (wallostonite, fibres courtes), dispersées dans la matrice. Les
composites à renforts particulaires sont peu utilisés à des fins structurelles, mais plutôt pour
des applications spécifiques ou des produits de grande consommation.
• les renforts à fibres longues : ces renforts se présentent sous la forme de fibres de faible
longueur (quelques dizaines de mm) orientées ou non, assemblées en nappes appelées mats.
Cet assemblage se fait généralement à l’aide d’un liant ou d’une couture.
• les renforts à fibres continues : ce type de renfort est constitué d’un assemblage ordonné de
fibres continues, c’est-à-dire dont la longueur est voisine des dimensions de la pièce finale.
Les directions privilégiées des fibres confèrent à la pièce finale ses caractéristiques
mécaniques principales, et le comportement de la matrice est de second ordre dans ces
directions. Le comportement de la matrice dans les autres directions, et la cohésion qu’elle
apporte, ont néanmoins une importance majeure.
Les fibres, longues ou continues, peuvent être formées à partir de nombreux matériaux : verre,
carbone, métaux, ou encore matériaux d’origine végétale (cellulose, lignine) ou animale (collagène).
Dans cette étude nous nous intéresserons exclusivement aux renforts à fibres continues. Dans ces
renforts l’assemblage des fibres peut être réalisé de différentes façons
• en couches de fibres unidirectionnelles (UD) superposées. Ces couches sont généralement
appelées « plis » et peuvent avoir des orientations de fibres différentes. Plusieurs plis
d’orientation différente peuvent être cousus : on obtient alors un NCF (Non Crimp Fabric)
• par tricotage
• par tressage
• par tissage avec un métier à tisser, en utilisant les armures classique de l’industrie textile ou
des armures plus complexes
Nous étudierons ici les renforts tissés : les fibres sont assemblées en paquets de fibres appelés
mèches ou torons, puis tissées selon des armures plus ou moins complexes. Dans les applications
structurelles, des renforts de fibres de verre et de carbone sont généralement utilisés. Le coût élevé de
ces matériaux les destine plutôt aux produits à haute valeur ajoutée pour lesquels l’allégement
représente un enjeu économique fort (e.g. pour l’industrie aéronautique/aérospatiale ou automobile).
Certaines applications « grand public » commencent néanmoins à voir le jour, souvent autant pour
l’esthétique de ces matériaux que pour leur caractéristiques mécaniques. Notons également que des
recherches sont menées à l’heure actuelle sur la possibilité de renforcer le renfort lui-même à l’aide de
nanotubes de carbone [GOD09, LOM11]. Ces nanotubes, destinés à des applications de haute
technicité, sont « cultivés » directement au sein du renfort, et permettent notamment d’améliorer le
comportement en rupture du composite.
Fibres de verre
Selon l’usage qui sera fait de la fibre, plusieurs types de compositions peuvent être utilisées comme
matière première. Parmi les différentes compositions classiques du verre, les plus couramment
utilisées pour former des fibres de verres destinées à des applications structurelles sont décrites dans le
Tableau 1.
VERRE FONDU
Four filière
Vitesse de refroidissement vT
(Température T, Diamètre d)
Organe d’ensimage
Organe(s) d’assemblage
Bobine
(Vitesse angulaire ω)
Figure 1-1. Procédé de filage et d'assemblage des mèches de fibres de verre [DAL08]
Fibres de carbone
Les fibres de carbone sont quant à elles obtenues à partir de tresses d’une matière première appelée
précurseur. Un procédé de carbonisation à haute température et sous atmosphère inerte transforme ces
tresses en fibres de carbone. Trois précurseurs peuvent être utilisés :
• les fibres de rayonne (cellulose régénérée),
• les fibres de brai obtenues par filage de résidu de houille ou de pétrole,
• le polyacronitrile (PAN).
Les fibres obtenues à l’aide de ces différents précurseurs sont respectivement appelées ex-rayonne,
ex-brai et ex-PAN. Parmi ces précurseurs, le premier n’est plus utilisé que pour l’obtention de fibres
de carbone à usage thermique. Si ce précurseur présente l’avantage de ne nécessiter aucun traitement
avant carbonisation, son faible rendement de carbonisation et le coût élevé du procédé de fabrication
associé l’ont rendu obsolète pour les applications structurelles.
Les fibres de brai permettent quant à elles d’obtenir des fibres de carbone possédant un très bon
taux de carbone, mais la difficulté d’élimination des impuretés dans le brai d’origine confère à ces
fibres une résistance à la traction relativement faible.
Enfin, les fibres de carbone peuvent être obtenues à partir de fibres acryliques, qui constituent le
précurseur le plus utilisé de nos jours. La figure 1-2 explique les différentes étapes de ce procédé.
Après oxydation et carbonisation des fibres de PAN, des fibres dites de haute résistance (HR) sont
obtenues. Si un procédé de graphitisation est mis en œuvre après la carbonisation (afin de parfaire la
structure moléculaire du carbone graphite), des fibres dites haut module (HM) sont obtenues, mais
dont la résistance est plus faible. Des variations de l’étape de graphitisation permettent d’obtenir
différentes gammes de fibres, depuis les fibres haute résistance/faible module jusqu’aux fibres faible
résistance/haut module.
Bobine de PAN
Oxydation
Carbonisation
Graphitisation
Traitement de surface
Ensimage
HR HM
Figure 1-2. Procédé de fabrication des fibres de carbone à partir d'un précurseur PAN [DUP08]
Fils de trame
Fils de chaîne
façon que les textiles d’habillement. Comme cela a été mentionné en introduction de cette section,
seules les armures 2D et 2.5D seront étudiées dans ce manuscrit.
Armures 2D
En ce qui concerne les tissus 2D il existe trois armures principales :
• la toile ou taffetas (plain weave), l’armure la plus simple : chaque fil de chaîne passe
alternativement au dessus et en dessous de chaque fil de trame (figure 1-4.a).
• le sergé nxm (nxm twill) : le fil de trame passe au dessus de n puis en dessous de m fils de
chaîne en décalant d'un fil à chaque passage (figure 1-4.b).
• le satin : les points de liage de la chaîne et de la trame sont disséminés de façon à atténuer
l’effet de diagonale présent sur le sergé (figure 1-4.c).
a. b. c.
Figure 1-4. Les différents motifs de tissage 2D. (a) Taffetas (b) Sergé 2x2 (c) Satin.
Filet de
Résine Moule Renfort Membrane drainage Joint
Ce procédé est moins onéreux et plus aisé à mettre en œuvre que le RTM, mais il ne permet pas de
contrôler finement l’épaisseur de la pièce finale, et l’état de surface du côté de la membrane est
généralement moins bon que celui du côté du moule.
a. b. c.
Figure 1-8. Schéma de principe du processus RTM. (a) Formage du renfort fibreux (b) injection
et polymérisation de la matrice (c) démoulage de la pièce finale
L’étape de mise en forme du renfort est une étape clé de ce procédé. La capacité du renfort à se
déformer pour prendre la forme du moule a en effet une influence majeure sur l’étape d’injection
(qualité de l’imprégnation, temps de remplissage du moule) ainsi que sur les caractéristiques
mécaniques de la pièce.
Le processus RTM classique est prévu pour fonctionner à basse pression (inférieure à 10 bar) et à
faible débit (de l’ordre du litre par minute). Certaines variantes ont néanmoins été développées pour
accélérer ou faciliter l’injection dans le cas de renforts peu perméables ou de résines très visqueuses :
une mise sous vide peut être utilisée (Vacuum Assisted RTM) ou encore une injection à haute pression
(High Speed RTM). Ces variantes rendent le procédé RTM adaptable à de nombreux types de renforts,
et donc à une gamme très large de pièces réalisables. Des pièces de grandes dimensions peuvent être
réalisées à l’aide de ce procédé.
Par ailleurs, le stockage de la matière première est généralement plus aisé pour le RTM que pour
les procédés nécessitant des préimprégnés (e.g. thermocompression), et la segmentation du procédé
permet un contrôle qualité à chacune des étapes. Ce procédé permet de réaliser de petites séries de
pièce à haute performance ou de grandes séries, s’adaptant ainsi aux exigences de l’industrie
aéronautique comme à celles de l’industrie automobile. Les moules des pièces complexes ou de
grande taille peuvent néanmoins nécessiter un effort de développement considérable. Ces nombreux
avantages incitent aujourd’hui les industriels à porter un grand intérêt au développement et à la
recherche associés à ce procédé. C’est dans ce cadre que se situe l’étude présentée dans ce manuscrit.
directement hérité de son comportement aux échelles inférieures. Trois échelles d’observation du
renfort sont généralement distinguées (figure 1-9) :
• échelle microscopique : échelle de la fibre
• échelle mésoscopique : échelle de la mèche
• échelle macroscopique : échelle de la pièce complète
À chacune de ces échelles le comportement mécanique ainsi que les défauts éventuels du renfort
sont différents. Cette section a pour but de décrire les mécanismes connus qui régissent le
comportement à chaque échelle, afin de pouvoir appréhender le comportement du renfort dans son
ensemble. Les moyens d’essai associés à ces mécanismes seront également présentés.
400 900
350 800
300 700
600
250
Effort (N )
Effort (N )
500
200
400
150 300
100 200
50 100
0 0
Si elle est clairement définie dans les matériaux continus, la notion de contrainte est moins évidente
pour une mèche de renfort tissé : considérons les sections apparentes S1 et S2 de la figure 1-11. Ces
deux sections contiennent le même nombre de fibres, seule la densité de fibres est différente.
Le nombre de fibres étant le même dans la section S1 que dans la section S2, le comportement de la
mèche en tension est identique dans les deux cas. C’est pourquoi définir une contrainte dans la mèche
comme l’effort divisé par la section apparente n’a que peu de sens, cette section n’ayant pas une
influence directe sur la réponse en effort du matériau. Le comportement en élongation de la mèche est
donc caractérisé non pas par un module mais par une rigidité, rapport de l’effort sur la déformation
dans la partie linéaire de la réponse en tension, et un effort maximum, tous deux exprimés en N.
F1 σ1
S1
S2
F2 σ2
Figure 1-11. Difficulté de la définition d'une contrainte, liée au changement de section apparente
de la mèche
rigidité de compaction de la mèche tend vers la rigidité de compaction du matériau constituant les
fibres. Ce comportement en compaction est difficile à caractériser expérimentalement, car la mise en
place d’essais de compaction pure sur une mèche sortie du tissu est difficilement envisageable. Dans
les études précédentes décrivant le comportement de la mèche à l’échelle mésoscopique, le
comportement en compaction de la mèche est généralement identifié par une méthode inverse à partir
d’un essai à l’échelle macroscopique [GAS00, HAG04, BAD08a, POT08]. C’est également ce qui
sera fait dans ce manuscrit.
a
.
a.
b
b.
Figure 1-12. Vue en coupe (par tomographie X) du plan transverse de la mèche (a) tissu au repos
(b) tissu sous tension équibiaxiale [BAD08a]
La compaction est l’un des modes de déformation privilégiés de la mèche, et elle intervient dans la
plupart des sollicitations auxquelles le tissu peut être soumis. Directement liée à cette compaction, la
densité de fibres au sein de la mèche a une influence considérable sur la perméabilité locale de la
mèche, et donc sur l’étape d’injection du procédé RTM. L’identification du comportement en
compaction est donc importante lorsque l’objectif est de calculer la perméabilité locale du renfort.
Chaque fibre hérite de son matériau constitutif et de sa géométrie une raideur en flexion. En
supposant le matériau constitutif de la fibre linéaire orthotrope, cette raideur peut être quantifiée
relativement facilement. Comme nous l’avons vu précédemment, une raideur en cisaillement
transverse peut être attribuée à la mèche du fait du frottement entre les fibres et de la présence
d’ensimage. Cette raideur est généralement trop élevée (i.e. la cohésion des fibres au sein de la mèche
est trop forte) pour que le cisaillement transverse se propage jusqu’aux extrémités de la mèche testée,
c’est pourquoi la mèche flambe latéralement (figure 1-15).
Le comportement de la mèche en flexion dépend donc du cisaillement transverse (i.e. de la
déformation locale) et de la courbure des fibres au sein de la mèche (i.e. de la courbure locale de la
mèche). Cette dépendance à la courbure locale permet d’affirmer que le comportement de la mèche en
flexion n’est pas un comportement de milieu continu classique (i.e. de Cauchy) : le comportement
local étant dépendant des rotations locales au sein du milieu, la mèche est en fait un milieu
micropolaire [COS09, MIN65, IND08]. Un modèle prenant en compte explicitement l’influence de la
courbure sur le comportement du renfort est proposé dans [MAG01].
La flexion de la maille élémentaire ne sera pas étudiée dans l’étude mésoscopique proposée dans ce
manuscrit, cette étude sera limitée aux chargements du renfort dans le plan. L’embuvage des mèches
étant très faible, la courbure de la mèche ne varie que très peu au cours de ces chargements, c’est
pourquoi cet aspect « de Cosserat » du comportement de la mèche ne sera pas pris en compte dans les
simulation à l’échelle mésoscopique.
Figure 1-15. Mécanismes intervenant dans la flexion d'une mèche de renfort tissé.
L’hypothèse de non glissement relatif des réseaux, qui sera faite dans le reste de ce manuscrit, est
donc vérifiée dans la plupart des cas de mise en forme classiques, mais il faut garder à l’esprit que
certains cas de chargement et certains tissus lui sont incompatibles.
Ces deux mécanismes induisent une forte non-linéarité du comportement du renfort en tension : sa
rigidité est faible tant que l’embuvage n’a pas été résorbé dans la direction de sollicitation, puis vaut la
somme des rigidités des n mèches de cette direction.
l0 ( 1+e) +d 140
120
l0 0
La forte réduction de l’embuvage dans la direction de sollicitation induit une forte augmentation de
l’embuvage dans la direction orthogonale. La longueur de la mèche restant constante, cette
augmentation de l’embuvage se traduit par un « avalement » des mèches dans cette direction. Ce
phénomène contredit l’hypothèse de non glissement entre les réseaux de chaîne et de trame formulée
au paragraphe [Link]. Cependant les renforts étudiés dans ce manuscrit ont pour la plupart un
embuvage si faible (inférieur à 1%) que ce phénomène d’avalement dans la direction orthogonale est
difficile à caractériser.
ε orth
k= (1.1)
ε obs
Un coefficient k = 1 signifie donc que les deux réseaux sont sollicités de façon identique (traction
équibiaxiale), et un coefficient k = 2 signifie que la déformation dans la direction du réseau observé
est deux fois plus faible que la déformation dans la direction de l’autre réseau.
uobs
Lobs
uobs
Le comportement en traction biaxiale est fortement non linéaire (figure 1-19), et cette non-linéarité
dépend, entre autres, du coefficient k : les mécanismes de déformation du tissu changent lorsque ce
coefficient varie. En traction équibiaxiale ( k = 1 ), c’est principalement la compaction de la mèche qui
accommode la déformation mais lorsque k = 0 ou k = +∞ , c’est la réduction de l’embuvage (et donc
le cisaillement transverse de la mèche) qui accommode la plus grande partie de la déformation.
Plusieurs dispositifs expérimentaux de traction biaxiale adaptés au tissu ont été conçus [KAW73a,
BUE01, WIL08], qui permettent de régler le coefficient de traction biaxiale (un exemple sur la figure
1-19). La fabrication des éprouvettes en croix utilisées dans ces dispositifs est un point clé de l’essai :
la raideur des mèches en tension étant très importante, et les déplacements très faibles, les mèches
doivent être parfaitement alignée dans les deux directions afin d’entrer en tension simultanément.
120
M èche seule
100
K=2
Effort par mèche (N)
80
K=1
60 K=0,5
Libre
40
20
Figure 1-19. Réseau de courbes de traction biaxiale pour un taffetas de verre [BUE98]
γ
2
l L
Figure 1-20. Cinématique et dispositif expérimental de l’essai de cadre (picture frame test).
π 1 d
γ= − 2 arccos + (1.2)
2 2 L
Afin que le comportement observé soit indépendant des dimensions de l’éprouvette, il faut
adimensionner l’effort F mesuré. Cela permettra également de comparer les résultats de cet essai
avec ceux d’autres essais de cisaillement. Une grandeur appelée couple surfacique C s permet une
mesure adimensionnée des efforts de cisaillement au sein du tissu [DEL09]. Cette grandeur permet de
comparer les essais réalisés avec différents tissus et/ou différents dispositifs d’essai. Elle est définie de
la façon suivante :
α
sin
L
Cs = F (1.3)
l2 2
La figure 1-21 montre une courbe obtenue à l’aide de cette approche pour un taffetas de verre.
L’essai de cadre, bien que théoriquement bien adapté à la caractérisation du cisaillement, présente
quelques difficultés :
• les mèches désalignées par rapport au côté du cadre vont être tendues lors de l’essai, créant
des efforts parasites [MIL07]. Les mèches étant très rigides, et le comportement en
cisaillement très souple, ces efforts parasites peuvent s’avérer plus importants que l’effort de
cisaillement, et ainsi perturber la mesure de façon significative,
• dans le modèle théorique, les mèches sont en liaison pivot avec le bord du cadre. En pratique
il n’est pas envisageable de réaliser une liaison pivot entre chaque mèche et le cadre. La
mèche flambe donc localement [DUM03a, DUM03b, WIL08, ZHU07], créant des
déformations parasites difficiles à évaluer,
• lorsque l’éprouvette est montée dans le cadre, une légère tension lui est appliquée. La valeur
de cette tension influe sur l’interaction entre les deux réseaux de mèches, et donc sur le
comportement en cisaillement [LAU08, LOM06]. Les efforts de glissement seront donc plus
importants au sein d’une éprouvette initialement plus tendue.
0.6
0.5
Couple surfacique (N/mm)
0.4
0.3
0.2
0.1
0 10 20 30 40 50 60
Angle de cisaillement (deg)
l’éprouvette : une zone non cisaillée (en bleu), une zone cisaillée d’un angle γ (en orange) et une zone
dite « demi cisaillée », cisaillée d’un angle γ / 2 (en vert). Avec les hypothèses précédentes, un angle
de cisaillement théorique peut être calculé en fonction de l’allongement de l’éprouvette :
π 1 d
γ= − 2arccos 1 + (1.4)
2 2 L − l
De même que pour l’essai de cadre, il faut définir une mesure adimensionnée des efforts de
cisaillement. La présence des zones non cisaillées et demi-cisaillées rend la relation entre l’effort
global sur l’éprouvette et l’effort local de cisaillement plus complexe que pour l’essai de cadre. Un
bilan de puissances permet d’exprimer, pour le bias extension test, le couple surfacique introduit
précédemment [DEL09] :
Cs (γ ) =
1 α 1 γ
F D sin − S 2C s (1.5)
S1 2 2 2
Dans l’expression (1.5), S1 et S2 sont respectivement l’aire initiale de la zone cisaillée et l’aire de
la zone demi cisaillée. Le premier terme est semblable à l’expression obtenue pour l’essai de cadre, et
le second est une correction liée à la présence des zones demi cisaillées. A partir de cette expression, le
couple surfacique peut être déterminé numériquement par une méthode séquentielle.
d
l
D α
Cet essai est exploitable tant que les hypothèses précédentes sont satisfaites. Dans les tissus peu
cohérents, c'est-à-dire dans lesquels les réseaux interagissent peu, il sera difficile d’appliquer cette
méthode car le glissement entre les deux réseaux apparaît dès le début de l’essai. A l’aide d’un outil de
corrélation d’images, il a notamment été montré que cet essai n’est plus en accord avec l’hypothèse de
non glissement des mèches au-delà d’un certain angle de cisaillement [DEL09].
0.6
PRESSIO N 1
0.5
PRESSIO N 2
Pression (M pa)
0.4
PRESSIO N 3
0.3
PRESSIO N 4
0.2 PRESSIO N 5
0.1
La figure 1-23 montre l’allure des courbes d’écrasement obtenues lors de ces essais : la raideur
initiale du renfort est très faible, puis elle augmente fortement lorsque la densité de fibres devient
importante. Plusieurs écrasements successifs de la même éprouvette ont été réalisés, illustrant ainsi
l’influence de l’histoire de la déformation sur le comportement du renfort. Bien que les essais
d’écrasement d’un renfort semblent assez simples à réaliser, plusieurs facteurs rendent
l’expérimentation difficile :
• le frottement du tissu avec les outils de compression ainsi que la taille de ces outils ont une
influence significative sur le comportement du renfort.
• le réarrangement des fibres lors de l’écrasement du renfort induit des phénomènes de
relaxation difficiles à quantifier, et qui rendent les résultats très sensibles à la vitesse de l’essai
et à l’histoire de la déformation. Ces phénomènes existent dans tous les essais mais leur
influence est particulièrement significative dans les essais d’écrasement.
Renfort
Plan
horizontal Plan
incliné
Si cet essai est très simple à réaliser, le dépouillement indiqué par la norme ASTM repose sur des
hypothèses incompatibles avec le comportement des renforts tissés : la raideur en flexion du tissu est
supposée constante. Or il a été montré que le comportement de flexion d’un renfort peut s’avérer
fortement non-linéaire [LIV64, KAW80, DEB10]. Le standard ASTM permet donc au mieux
l’obtention d’un ordre de grandeur de la rigidité en flexion du renfort, ne prenant pas en compte cette
non linéarité.
Le flexomètre de De Bilbao
Ce dispositif ressemble au dispositif décrit par la norme ASTM : le tissu repose sur des lames qui
sont retirées progressivement (figure 1-26), entraînant une flexion du tissu sous l’action de son propre
poids. Un système de mesure optique permet, à l’aide d’un programme de traitement d’images,
d’obtenir le profil de l’échantillon pour chaque longueur de surplomb. Un dépouillement complexe
basé sur l’ensemble des profils de flexion obtenus est alors mis en place, qui permet de prendre en
compte la non-linéarité du comportement en flexion du renfort. Ce dispositif, et la méthode de
dépouillement associée, permettent une analyse plus fine du comportement du renfort que le standard
ASTM, et est mieux adapté aux renforts de composites que le KeS-FB2 car il permet de mesurer des
raideurs de flexion plus importantes.
Dans la littérature, de nombreux articles portent sur la simulation des cordes et des câbles, qui
s’apparentent à des mèches comportant moins de fibres et fortement torsadées [GHO07, JIA99,
NAW00, SRE06]. Les études spécifiques aux renforts tissés sont plus rares. Une modélisation du
tissage a notamment été développée [ZHO04], qui permet de prévoir, dans certaines limites, la
géométrie des mèches au sein de la maille élémentaire du renfort. Des modèles de mailles élémentaires
complètes ont été développés [DUR10], qui permettent notamment de visualiser, au cours d’une
sollicitation, le réarrangement des fibres au sein de la mèche, et les contraintes dans chaque fibre (voir
figure 1-27).
1
Dans ce manuscrit, le terme plissement sera préféré au terme pli afin d’éviter la confusion avec le terme
« pli » désignant généralement une couches au sein d’un renfort multicouches.
La limitation principale de ce type de modèles est le temps de calcul nécessaire : à l’heure actuelle,
seules quelques dizaines de fibres par mèche peuvent être modélisées, alors qu’une mèche de renfort
tissé peut en contenir quelques dizaines de milliers. Ceci est très limitant car lorsque les fibres sont
rassemblées par paquets, car dans ce cas le comportement de ces paquets est plus proche de celui de la
mèche, beaucoup plus complexe que de celui des fibres.
Ces hypothèses éloignent toute considération mécanique, et donc toute différentiation entre les
différents renforts : le résultat de l’algorithme du filet sera identique quel que soit l’armure du renfort
considéré ou la raideur de ses mèches. Par ailleurs, l’absence de frottements avec les outils de mise en
forme empêche la prise en compte d’éventuels serre-flans.
Dans cette approche, le renfort est donc semblable à un ensemble de barres articulées non
déformables (figure 1-28). Partant d’une ligne initiale, l’algorithme calcule la position des points
connexes sur des géodésiques de la surface à draper issues de ce point. Dans l’étude des renforts tissés,
l’algorithme du filet peut être utilisée comme approche préliminaire, permettant de déterminer le degré
de complexité du formage d’une pièce.
Figure 1-29. Différents types de raideurs des modélisations mésoscopiques discrètes [SZE05]
a. a
.
b
b.
.
Figure 1-30. Simulation avec un modèle hypoélastique (a) et expérience (b) du drapage d'un
double dôme [KHA10]
A ce jour, plusieurs lois de comportement ont été développées à l’aide de différentes approches,
notamment des approches hypoélastiques [PEN05, YU05, KHA10] et des approches hyperélastiques
[TEN07, AIM10, DRI10]. Un exemple de résultat obtenu à l’aide de l’une de ces méthodes est
présenté sur la figure 1-30.
Enfin, les approches semi discrètes sont une combinaison des deux approches précédentes : des
éléments finis sont formulés dans lesquels la puissance des efforts internes est subdivisée en
différentes contributions, liées aux différentes rigidité du tissu. Deux modèles de ce type ont été
développés, le premier dédié aux tissus 2D utilisant des éléments de coque ou de membrane [BOI06,
HAM07b, HAM09], et le second dédié aux interlocks 3D utilisant des éléments 3D [DEL09]. Des
résultats d’emboutissages de renfort obtenus par ces deux modèles sont présentés sur la figure 1-31.
a. b.
Figure 1-31. Schématisation des éléments formulés et résultats de simulation de mises en formes
obtenus pour le modèle semi discret 2D (a) et le modèle semi discret 3D (b)
Une approche analytique intéressante a également été développée par Lomov [LOM06]. Cette
approche est basée sur une minimisation de l’énergie de déformation des mèches au sein du renfort.
Un système d’équations est construit à partir d’une paramétrisation de la géométrie des mèches, et de
différentes contributions mécaniques : tension des mèches, flexion des mèches, torsion des mèches,
compression des mèches, contacts aux intersections de mèches, contacts latéraux entre mèches et
frottement entre mèches. De petites portions de mèches (elementary bent intervals) sont définies pour
discrétiser le problème de minimisation, qui est alors résolu afin d’obtenir l’état d’équilibre mécanique
du renfort. Une approche similaire basée sur un principe de stationnarité de l’énergie potentielle totale
du système traduisant l’équilibre entre les efforts extérieurs appliqués et la réponse du système a
également été proposée [SAG03].
Les approches analytiques fonctionnent bien pour déterminer le comportement mécanique de la
maille élémentaire dans en tension uniaxiale et en tension biaxiale. L’application de ces approches à
des cas plus complexes tels que le cisaillement est possible mais difficile, car la simplification de la
géométrie de la mèche rend les conditions de contact latéral entre les mèches peu réalistes. Face aux
limitations de ces approches analytiques, les simulations par éléments finis apparaissent comme une
très bonne alternative. Les simulations par éléments finis à l’échelle mésoscopiques nécessitent trois
données principales : la description de la géométrie de la maille élémentaire par des éléments finis
(3D), la loi de comportement de la mèche, et les conditions aux limites appliquées à la maille
élémentaire. Ces différents aspects sont complexes et ont fait l’objet de plusieurs études, qui sont
présentées dans cette section.
• l’utilisation d’un tomographe à rayons X : cette méthode, plus récente, nécessite un matériel
sophistiqué mais ne présente pas les défauts des deux méthodes précédentes [BAD08b].
Par ailleurs, le modèle géométrique obtenu doit être cohérent : il doit garantir l’absence
d’interpénétration entre les mèches. Cette condition est particulièrement délicate aux points de
croisement des mèches. La difficulté principale de la modélisation géométrique du renfort réside donc
dans la nécessité de décrire la géométrie de la mèche avec une complexité suffisante, tout en assurant
la cohérence du modèle.
Plusieurs modèles ont été proposés dans la littérature pour déterminer cette géométrie de la maille
élémentaire. Une premier, mis en œuvre dans le logiciel Wisetex, est développé à l’université de
Louvain [LOM00, VER05, LOM06]. Il consiste à définir une paramétrisation de la géométrie de la
maille élémentaire, puis à déterminer ces paramètres en minimisant l’énergie de flexion des mèches de
la maille élémentaire. La géométrie est décrite par la ligne moyenne des mèches, décomposée en
segments ondulés élémentaires (paramétrisés par des polynômes d’ordre 5), et par les sections des
mèches à chaque croisement, de formes lenticulaires ou elliptiques (paramétrisées par deux
dimensions caractéristiques). Un effort transversal est défini à chaque croisement, permettant de
décrire l’interaction entre les deux réseaux de mèches. Cet effort transversal est relié au comportement
de la mèche en compaction, et ainsi à la forme des sections transverses. Le problème de minimisation
de l’énergie est alors résolu afin de déterminer l’ensemble des paramètres de la description
géométrique de la maille élémentaire. Ce modèle permet de déterminer rapidement des géométries de
renforts très complexes, y compris des renforts 2,5D et 3D. Il permet également de prendre en compte
les ondulations latérales des mèches qui apparaissent notamment lorsque la géométrie des mèches
n’est pas symétrique au niveau des points de croisement, ce qui est le cas dans la plupart des armures
autres que le taffetas. Il s’agit également de l’un des seuls modèles capables de fournir des
informations sur l’interaction entre les deux réseaux de mèches ainsi que l’orientation des fibres dans
un tissu au repos. Initialement ce modèle ne garantissait pas la condition de cohérence : des
intersections entre les mèches étaient susceptibles d’apparaître. Une méthode a été implémentée par la
suite afin résoudre ces problèmes d’incohérence [VER05]. Une transition directe entre les modèles
obtenus dans Wisetex et des calculs par éléments finis de comportement mécanique ou de perméabilité
de la maille élémentaire de renfort est maintenant possible.
Une alternative à Wisetex est le logiciel Texgen, développé à l’université de Nottingham [SHE07].
A l’inverse de Wisetex, Texgen ne fait pas intervenir de considérations mécaniques dans le calcul de la
géométrie du renfort au repos. Texgen ne permet pas non plus de calculer des géométries déformées de
maille élémentaire. Il permet toutefois de réaliser facilement des modèles de maille élémentaires
complexes (2D, 2.5D, 3D, NCD, tresses, tricots), dans lesquels les mèches peuvent avoir des sections
de formes complexes et variables le long de la ligne moyenne. La gestion de l’interpénétration entre
les mèches est ici aussi une question majeure. Notons qu’avec Texgen il est possible d’exporter
directement le maillage de la maille élémentaire dans un format de mise en donnée pour Abaqus.
Le dernier type de modèle qui a été envisagé pour la modélisation de la géométrie de la maille
élémentaire est un modèle géométrique proposé par Gilles Hivet [HIV05, HIV08]. Ce modèle ne
s’intéresse qu’aux renforts tissés 2D mais propose une modélisation géométrique garantissant
l’absence d’interpénétrations entre les mèches, c’est pourquoi il sera utilisé dans ce manuscrit pour
générer les modèles de mailles élémentaires. Dans ce modèle, la section de la mèche peut varier le
long de la ligne moyenne. Les sections sont définies en un certain nombre de points de contrôle sur
cette ligne moyenne, nombre variable selon le type d’armure considéré (figure 1-33). Ces points de
contrôle définissent les différentes zones caractéristiques de la mèche : les zones sans contact avec
d’autres mèches (e.g. M3-M4), les zones avec contact courbe (e.g. M2-M3) et les zones avec un
contact plan (e.g. M1-M2).
Figure 1-33. Définition des points de contrôle de la ligne moyenne d'un sergé 3x2 [HIV05]
Les sections de la mèche sont définies en chacun de ces points de contrôle, c'est-à-dire pour
chacune des différentes zones. Le contour de chaque section est divisé en quatre parties, modélisées
par quatre coniques (figure 1-34). Une dégénérescence des coniques situées sur les bords de la mèche
permet d’obtenir des formes simplifiées dont le bord est réduit à un point ou à un segment, ce qui
pourra par la suite faciliter la discrétisation de la géométrie en éléments finis. Notons que pour les
armures taffetas les sections sont symétriques (figure 1-34b).
La géométrie de la mèche entre les sections est obtenue par un balayage le long de la ligne
moyenne, avec interpolation des sections des points de contrôle. En utilisant les formes simplifiées, il
a été montré que trois paramètres sont nécessaires pour modéliser un taffetas équilibré, et six pour un
taffetas déséquilibré. Un exemple de maillage d’un taffetas de verre issu de cette modélisation
géométrique est présenté sur la figure 1-35.
alignées et de vides entre ces fibres. Ce comportement doit rendre compte des spécificités d’un tel
assemblage, liées notamment de la très forte anisotropie de sa structure et aux possibles glissements
entre les fibres au sein de la mèche. La définition de cette loi de comportement n’est donc pas triviale.
La plupart des études précédentes utilisent des modèles hypoélastiques pour décrire le
comportement de la mèche. L’un des premiers modèles hypoélastiques pour la mèche a été proposé
par Gasser dans [GAS00] : la mèche est considérée comme un matériau orthotrope, et seules les
contributions de la tension de la mèche dans la direction des fibres, et de sa compaction dans le plan
transverse sont pris en compte.
Par la suite, les travaux de Hagège [HAG04] ont montré que l’utilisation des dérivées classiques de
Jaumann et de Green-Naghdi ne permet pas de garantir l’objectivité des lois de comportement
hypoélastique pour la mèche : les repères associés à ces dérivées tournent avec des rotations moyennes
de la matière, et donc ne suivent pas exactement la direction matérielle liée aux fibres. Les travaux de
Hagège introduisent donc la définition d’une nouvelle dérivée, objective, calculée dans un repère
suivant exactement la direction matérielle des fibres. Ces travaux ont mené notamment à une
description satisfaisante du comportement en traction biaxiale d’un renfort tricoté.
Ces travaux ont été poursuivis par Badel [BAD08a] qui a appliqué cette loi de comportement
orthotrope à des renforts tissés, et qui a notamment mis au point une simulation satisfaisante du
cisaillement d’une maille de renfort tissé (figure 1-36).
Ces lois de comportement hypoélastiques posent néanmoins des problèmes d’identification, liés à
la difficulté de réaliser des essais à l’échelle de la mèche. Par ailleurs, si le formalisme de la dérivée
objective suivant la direction des fibres apparaît aujourd’hui relativement simple, l’extension de ce
formalisme au suivi simultané de plusieurs directions de fibres, initié dans [HAG04], est relativement
ardue. Un tel formalisme serait pourtant intéressant pour étendre les modèles mésoscopiques aux tissus
interlocks 2,5D.
Ces modèles hypoélastiques posent également un problème : ils ne sont pas élastiques, et
l’anélasticité introduite est difficile à maîtriser car elle dépend du code de calcul utilisé. Ainsi, il faut
s’attendre à dissiper de l’énergie lors d’un chargement cyclique sur la mèche, bien qu’aucun terme de
dissipation n’ait été introduit dans la loi de comportement.
Conclusion du chapitre 1
Dans ce premier chapitre, les renforts tissés de composites ont été présentés dans leur globalité, de
la fibre et sa fabrication au renfort final avec son armure. Cette description de la structure du renfort
aux différentes échelles nous a permis de mieux comprendre les mécanismes de déformation du
renfort. Le procédé RTM a également été présenté, ainsi que les enjeux industriels associés à la
simulation de ce procédé.
Les essais « classiques » utilisés pour décrire le comportement mécanique du renfort ont été
présentés : tension uniaxiale et biaxiale, picture frame test, bias extension test, flexion, écrasement.
Ces essais seront utilisés dans la suite de ce manuscrit pour identifier ou valider les lois de
comportement proposées.
Finalement, un état de l’art de la simulation de la mise en forme des renforts tissés a été réalisé,
présentant les modèles micro-, méso- et macroscopiques existant dans la littérature, et permettant de
simuler la déformation du renfort aux différentes échelles.
A l’échelle mésoscopique, les lois de comportement hypoélastiques existantes sont difficiles à
identifier et ne permettent pas une analyse fine des énergies mises en jeu eu sein de la maille
élémentaire du renfort. Par ailleurs, le formalisme lié au suivi de la direction des fibres, utilisant des
dérivées objectives, est assez complexe, c’est pourquoi des hypothèses sont généralement réalisées qui
ne permettent pas d’assurer un suivi exact de la matière. Une nouvelle approche sera présentée dans ce
manuscrit pour décrire le comportement de la mèche de renfort tissé, utilisant une loi de comportement
hyperélastique.
Aucune approche continue 3D n’existe dans la littérature pour simuler le comportement des
interlocks épais. Nous verrons que le formalisme hyperélastique permettra d’étendre la loi de
comportement de la mèche à une loi de comportement adaptée à des matériaux à plusieurs directions
de fibres, permettant ainsi de simuler la mise en forme des renforts interlocks.
Afin de formuler ces nouvelles lois de comportement, un aperçu du formalisme de la mécanique
des grandes transformations et de l’hyperélasticité anisotrope sera proposé au chapitre suivant. La
méthode qui sera utilisé pour implémenter ces lois de comportement dans le code de calcul par
éléments finis Abaqus/Explicit sera également présentée.
Dans ce chapitre, les concepts et les méthodes nécessaires à la mise en place de lois de
comportement hyperélastiques anisotropes non-linéaires en grandes transformations et leur
intégration dans le code Abaqus/Explicit seront présentés. Les concepts de base de la mécanique des
milieux continus en grandes transformations seront brièvement rappelés. Le formalisme de
l’hyperélasticité isotrope et anisotrope pourra alors être introduit, et les modèles hyperélastiques les
plus classiques seront présentés succinctement. Les différents types de schémas d’intégration
temporelle seront également décrits, et le choix du code de calcul Abaqus/Explicit pour la simulation
de la mise en forme des renforts tissés de composites sera justifié. La méthodologie d’implémentation
d’une loi de comportement dans ce code de calcul, sous la forme d’une subroutine utilisateur VUMAT,
sera présentée. Nous aurons alors à notre disposition tous les outils nécessaires au développement et
à l’implémentation de lois de comportement hyperélastiques adaptées au comportement mécanique
des renforts tissés de composites.
Sommaire
x = χ (X ) (2.1)
dx = xq − x p = χ ( X P + d X ) − χ ( X P ) (2.2)
∂χ
dx = ⋅d X = F ⋅d X (2.3)
∂X
Ainsi défini, le tenseur gradient de la transformation F constitue une description locale (au
premier ordre) de la transformation du solide.
L’application χ étant bijective, son jacobien J sera toujours strictement positif, ce qui revient à
dire que le volume de l’élément de matière défini par le volume élémentaire dV ne peut s’annuler.
2
Dans toute cette partie, et sauf mention contraire, les grandeurs écrites en minuscule seront relatives à la
configuration actuelle tandis que les grandeurs écrites en majuscule seront relatives à la configuration initiale.
dan = J F −T ⋅ N dA (2.5)
F = R ⋅U = V ⋅ R (2.6)
où R est un tenseur de rotation (i.e. orthogonal direct) et où U et V sont des tenseurs symétriques
définis positifs appelés tenseurs de déformation pure. Afin d’étudier les variations de longueur et
d’angle au sein du solide, on peut étudier le produit scalaire de deux vecteurs élémentaires en
configuration déformée :
= dU ⋅ F T ⋅ F ⋅ dV (2.7)
= dU ⋅ C ⋅ dV
du ⋅ dv − dU ⋅ dV = dU ⋅ C ⋅ dV − dU ⋅ dV
( )
= dU ⋅ C − I ⋅ dV (2.8)
= dU ⋅ 2 E ⋅ dV
( )
Le tenseur E = 1/ 2 C − I introduit ainsi est appelé tenseur des dilatations de Green-Lagrange :
en petites déformations, ce tenseur se confond avec le tenseur des déformations linéarisé ε . Les deux
tenseurs précédents sont définis sur la configuration initiale. A ce titre, il sont généralement dits
« lagrangiens ». On peut définir de la même façon leurs équivalents en configuration actuelle, qui
seront dits « eulériens »3 :
b = F ⋅ FT =V 2 (2.9)
3
Un tenseur sera dit lagrangien (resp. eulérien) lorsque ses quantités d’entrée et de sortie sont exprimées
dans la configuration initiale (resp. actuelle). Un tenseur dont la quantité d’entrée et la quantité de sortie ne sont
pas exprimées dans la même configuration sera dit bipoint.
a=
1
2
(
I − b−1 ) (2.10)
Les tenseurs définis par (2.9) et (2.10) seront appelés respectivement tenseur de Cauchy-Green
gauche et tenseur d’Euler-Almansi.
d xɺ = Fɺ ⋅ d X = Fɺ ⋅ F −1 ⋅ d x = l ⋅ d x (2.11)
Le tenseur l est appelé tenseur gradient de vitesse. Il est classiquement décomposé en un terme
symétrique d (tenseur taux de déformation) et un terme antisymétrique ω (tenseur taux de rotation) :
l = d +ω
(2.12)
d=
1
2
( )
1
F + FT , ω = F − FT
2
( )
Le tenseur taux de déformation est une mesure de la vitesse de déformation dans la configuration
actuelle. En différentiant l’équation (2.7), on peut également obtenir une mesure de la vitesse de
déformation dans la configuration initiale :
1
Eɺ = Cɺ = F T ⋅ d ⋅ F (2.13)
2
df
t= (2.14)
ds
t = σ ⋅n (2.15)
dS
ds
N
dF
n
F df
d f = σ ⋅ n ds (2.16)
Le tenseur de Cauchy représente donc les efforts internes exprimés dans la configuration actuelle.
De même que pour les différentes mesures de déformation établies précédemment, il est possible de
définir d’autres mesures des efforts internes. Ainsi, en notant d F , dS et N l’effort, la surface et sa
normale ramenés en configuration initiale, on définit les trois tenseurs S , P et τ :
d f = P ⋅ N dS
d F = S ⋅ N dS (2.17)
d F = τ ⋅ n ds
Ces tenseurs sont respectivement appelés premier tenseur de Piola-Kirchhoff, second tenseur de
Piola-Kirchhoff (ou tenseur de Boussinesq) et tenseur de Kirchhoff. En petites déformations, ces
quatre tenseurs sont identiques. Les équations (2.3) et (2.5) permettent de démontrer les relations
suivantes :
τ = Jσ = F ⋅ S ⋅ F T = P ⋅ F T (2.18)
Avec la définition précédente des adjectifs lagrangien et eulérien, le tenseur des contraintes de
Cauchy σ , qui permet d’obtenir le vecteur contrainte dans la configuration actuelle en fonction de la
normale dans cette même configuration, sera dit eulérien. A l’inverse, le second tenseur des
contraintes de Piola-Kirchhoff S , qui permet d’obtenir le vecteur contrainte dans la configuration
initiale en fonction de la normale dans cette même configuration, sera dit lagrangien. Le premier
tenseur de Piola-Kirchhoff P et le tenseur de Kirchhoff τ ne sont ni lagrangiens ni eulérien : ils sont
appelés tenseurs bipoints.
Ω
(
Pi = − ∫ σ : grad ( v ) ) T
dΩ (2.19)
∂v i ∂v ∂X j ∂ ∂x i ∂X j
grad ( v ) = = i = = Fɺ ⋅ F −1 = l (2.20)
∂x j ∂X k ∂x j ∂t ∂X k ∂x j
Pi = − ∫ σ : d d Ω (2.21)
Ω
Pi = − ∫ ρ pi d Ω = − ∫ ρ0 pi d Ω0 (2.22)
Ω Ω0
En utilisant l’expression (2.21) de la puissance des efforts intérieurs, ainsi que des relations (2.18)
entre les différentes mesures de contraintes, cette puissance massique s’écrit :
1 1 1
pi = σ :d = S : Eɺ = P : Fɺ (2.23)
ρ ρ0 ρ0
ce qui montre la dualité entre les autres mesures de déformation et d’efforts internes introduites
précédemment. Le terme σ : d représente la puissance des efforts internes par unité de volume final,
tandis que les termes S : Eɺ et P : Fɺ représentent la puissance des efforts internes par unité de volume
initial.
( 1
Φ 0 = − ρ0 ψɺ − θɺη − Q ⋅
θ
∂θ
∂x
)+ S : Eɺ ≥ 0 (2.24)
( 1
Φ = − ρ ψɺ − θɺη − q ⋅
θ ∂X )
∂θ
+σ : D ≥ 0 (2.25)
avec Φ 0 et Φ les dissipations par unité de volume (initial et final), ψ l’énergie libre spécifique
(i.e. par unité de masse), θ la température et Q et q les vecteurs flux de chaleur en configuration
initiale et actuelle. En théorie « purement mécanique », il est courant de s’affranchir des effets
thermiques, c'est-à-dire de considérer la température homogène et la transformation isotherme. Les
dissipations s’écrivent alors :
Φ 0 = − ρ0ψɺ + S : Eɺ (2.26)
Φ = − ρψɺ + σ : D (2.27)
Un matériau hyperélastique est par définition un matériau non dissipatif (i.e. pour lequel les
dissipations (2.26) et (2.27) sont nulles) dont l’énergie de déformation par unité de volume initial
w = ρ0ψ ne dépend que de l’état actuel des déformations du matériau. Ces deux considérations se
traduisent par les équations suivantes :
wɺ ( F ) = S : Eɺ
(2.28)
w = w ( F ) = ρ0ψ
On peut montrer par ailleurs [CIA88] qu’afin que le matériau respecte le principe d’indifférence
matérielle, il faut et il suffit que cette énergie de déformation soit une fonction du tenseur de Cauchy-
Green droit. Nous écrirons donc directement les relations (2.28) sous la forme :
( )
w = w C = ρ 0ψ
(2.29)
( ) 1
wɺ C = S : Eɺ = S : Cɺ
2
∂w 1 ∂w
( )1
wɺ C = S : Cɺ ⇔
2
− S : Cɺ = 0 ⇔ S = 2
∂C 2 ∂C
(2.30)
Cette dernière équation constitue l’expression générale des lois de comportement hyperélastique.
On peut également la mettre sous la forme :
2 ∂w T
σ= F⋅ ⋅F (2.31)
J ∂C
Pour établir la loi de comportement d’un matériau hyperélastique, il suffit donc de déterminer
l’énergie de déformation w correspondant au comportement mécanique du matériau considéré. Afin
d’être admissible, cette énergie de déformation doit vérifier plusieurs conditions :
• w s’annule lorsque le matériau est dans son état naturel (i.e. sans sollicitation ni contraintes
résiduelles) :
( )
w I =0 (2.32)
( ) ( )
w F = w QF , ∀Q ∈ SO3 (2.33)
( ) ( )
w F = w FQ , ∀Q ∈ G ⊂ SO3 (2.34)
où G est le groupe de symétrie du matériau. La condition (2.32) dépend de la forme choisie pour
l’énergie de déformation, elle devra donc être vérifiée lors de l’établissement concret de la loi de
comportement. La condition (2.33) est satisfaite par l’utilisation du tenseur de Cauchy-Green droit
comme argument de w [CIA88]. En ce qui concerne la condition (2.34), nous verrons par la suite, à
travers les cas de l’isotropie et le l’isotropie transverse, que le respect des symétries du matériau est
satisfait par l’introduction d’invariants du tenseur de Cauchy-Green droit.
( ) ( )
w C est isotrope ⇔ ∃wɶ : ℝ 3 → ℝ, w C = wɶ ( I 1 , I 2 , I 3 ) (2.35)
Les invariants classiquement utilisés pour définir la fonction wɶ sont les coefficients du polynôme
caractéristique de C (ou les valeurs propres λ1 , λ2 , λ3 de F qui y sont liées directement):
( )
I 1 = T r C = λ12 + λ22 + λ32
1
( ( ) − Tr (C )) = ( λ λ ) + ( λ1λ3 ) + ( λ2 λ3 )
2
I2 = 2 2 2 2
Tr C 1 2 (2.36)
2
= Det (C ) = ( λ λ λ )
2
I3 1 2 3
Dans un souci de simplicité, la même notation sera utilisée pour w et wɶ . Différentes lois de
comportement isotropes peuvent être trouvées dans la littérature : ces lois sont classiquement utilisées
pour étudier le comportement de matériaux caoutchoutiques, ainsi que de certains tissus biologiques.
Lorsque le matériau étudié est considéré comme incompressibles (ce qui est souvent le cas dans les
matériaux caoutchoutiques notamment), une liaison interne est souvent postulée sous la forme I 3 = 1 ,
qui se traduit dans la loi de comportement par une pression hydrostatique p :
2 ∂w T
σ= F⋅ ⋅ F − pI (2.37)
J ∂C
∞
w= ∑ C (I − 3) ( I 2 − 3)
i j
ij 1 (2.38)
i , j =0
avec C 00 = 0 afin que l’énergie de déformation s’annule lorsque le matériau n’est pas sollicité.
Trois modèles classiques découlent de ce modèle de Rivlin :
• le modèle néo-hookéen : w = C10 ( I 1 − 3 )
• le modèle de Mooney-Rivlin : w = C10 ( I 1 − 3 ) + C 01 ( I 2 − 3 )
• le modèle de Yeoh : w = C10 ( I 1 − 3 ) + C 20 ( I 1 − 3 ) + C 30 ( I 1 − 3 )
2 3
Afin que ces modèles assurent une bijection entre les contraintes et les déformations, les
coefficients Cij doivent être positifs.
N
µk α
w =∑ ( λ1 + λ2α + λ3α − 3)
k k k
(2.39)
k =1 α k
On peut montrer que les modèles néo-hookéen et de Mooney-Rivlin sont également des cas
particuliers du modèle d’Ogden. Ce modèle est particulièrement adapté aux très grandes déformations
de matériaux caoutchoutiques.
1
F = J3F (2.40)
Les modèles précédents sont alors utilisés en utilisant les invariants (dits « réduits ») I 1 et I 2 de
F , et un terme en I 3 = J 2 permet d’ajouter au modèle l’énergie de déformation liée à la dilatation du
matériau.
( ) (
w C = wˆ C , M 1 ...M N ) (2.41)
Comme pour le cas isotrope, les théorèmes de représentation des fonctions isotropes à valeurs
réelles et à arguments scalaires, vectoriels et/ou tensoriels, donnent alors l’existence d’un nombre fini
de scalaires ( I 1 ,..., I n ) ( n dépendant du nombre et du type d’arguments de la fonction considérée)
⌢
ainsi qu’une fonction w tels que
( )
wˆ C , M i = w ( I 1 ,..., I n )
⌢
(2.42)
w it = w it ( I 1 , I 2 , I 3 , I 4 , I 5 ) (2.43)
∂w ∂w ∂I 1 ∂w ∂I 2 ∂w ∂I 3 ∂w ∂I 4 ∂w ∂I 5
S=2 = 2 + + + + (2.44)
∂C ∂I 1 ∂C ∂I 2 ∂C ∂I 3 ∂C ∂I 4 ∂C ∂I 5 ∂C
L’expression (2.44) est la forme la plus générale des lois de comportement hyperélastique pour les
matériaux isotropes transverse.
Un certain nombre de modèles isotropes transverses peuvent être trouvés dans la littérature. La
plupart de ces modèles décrivent le comportement de matériaux dans lequel une matrice et un renfort
peuvent être distingués : matériaux caoutchoutiques renforcés [AGO09, DIA04], composites
fibres/matrice organique [GUO07], tissus organiques comme les parois artérielles [BAS09, HOL00]
ou les ligaments [HIR00]. Une superposition du comportement de la matrice et du comportement du
renfort est généralement utilisée, en additionnant deux densités d’énergie distinctes pour chacune de
ces phases. Par ailleurs, ces modèles sont généralement définis en fonction des invariants classiques ou
des élongations du matériau, ce qui ne sera pas le cas pour l’approche qui sera proposée dans ce
manuscrit.
trame du renfort, et que le vecteur M 3 est initialement orthogonal à ces deux directions en
configuration initiale.
Le théorème de représentation permet alors la représentation de la fonction densité d’énergie de la
façon suivante [QUA94, ITS07] :
M uɺɺ + C uɺ = f ext
− f int
(2.46)
où u est une matrice colonne contenant tous les degrés de liberté du modèle (déplacements des
nœuds), M est la matrice de masse, C est la matrice d’amortissement, et f int et f ext sont des
matrices colonnes contenant respectivement les efforts internes et externes aux nœuds.
La résolution de cette équation dans le temps se fait à l’aide d’un schéma d’intégration temporelle.
Le temps est d’abord discrétisé en incréments espacés par un pas de temps, ce pas de temps pouvant
être constant ou variable selon le type de calcul effectué. L’équation (2.46) discrétisée en temps
s’écrit, à l’incrément i :
ext , ( i ) int, ( i )
M uɺɺ(i ) + C uɺ (i ) = f − f (2.47)
u ( i + 1) (i )
= u + ∆t uɺ (i ) (i )
+
( ∆t ) (i ) 2
uɺɺ(i ) + β ( ∆t (i ) ) uɺɺɺ (i )
2
2 (2.48)
uɺ (i +1) = uɺ (i ) + ∆t (i )uɺɺ(i ) + γ ( ∆t (i ) 2
) uɺɺɺ (i )
Il est ensuite supposé que l’accélération est linéaire sur l’incrément, ce qui permet d’écrire le jerk
sous la forme :
( i + 1)
uɺɺ − uɺɺ(i )
uɺɺɺ (i ) = (2.49)
∆t (i )
En remplaçant cette expression du jerk dans les formules (2.48), les expressions classiques du
déplacement et de la vitesse associées au schéma de Newmark sont obtenues :
( ∆t ) (i ) 2
u ( i + 1) (i )
= u + ∆t uɺ (i ) ( i )
+
2
( (1 − 2β ) uɺɺ (i )
+ 2 β uɺɺ(i +1) ) (2.50)
uɺ ( i + 1) (i )
= uɺ + ∆t (i )
( ( 1 − γ ) uɺɺ(i ) + γ uɺɺ( ) ) i +1
Les paramètres β et γ permettent de décliner les équations (2.50) et d’obtenir différents schémas
d’intégration, adaptables au problème donné. La stabilité et la précision du schéma dépendent
fortement de ces deux paramètres.
Lorsqu’il est possible de connaître le déplacement à l’incrément i + 1 en fonction du déplacement,
de la vitesse et de l’accélération à l’incrément i , le schéma est dit explicite. Par exemple, le couple
β = 0, γ = 0 donne un schéma purement explicite et le couple β = 0, γ = 1/ 2 donne un schéma
explicite aux différences centrées. A l’inverse, lorsque le déplacement à l’incrément i + 1 en fonction
du déplacement ne peut pas être calculé explicitement à partir de la vitesse et de l’accélération à
l’incrément i , le schéma sera dit schéma implicite. Une description succincte de ces deux types de
schémas est proposée ci-après.
Lorsque 2 β ≥ γ ≥ 1/ 2 , ces schémas sont inconditionnellement stables. Les pas de temps peuvent
donc prendre des valeurs assez grandes, permettant ainsi d’effectuer un calcul en peu d’itérations.
Dans le cas de matériaux ou de structures au comportement non-linéaire, la convergence de
l’algorithme permettant de satisfaire l’équation de la dynamique à chaque incrément peut néanmoins
s’avérer difficile, auquel cas une diminution du pas de temps est nécessaire, augmentant d’autant le
temps de calcul. Pour des structures au comportement fortement linéaire (matériaux avec un
raidissement exponentiel, conditions de contact sévères, …) une forte diminution du pas de temps peut
s’avérer nécessaire, pouvant mener à des temps de calcul prohibitifs.
Schéma β γ
Accélération moyenne 1/ 4 1/ 2
Accélération linéaire 1/ 6 1/ 2
Fox-G oodwin 1/12 1/ 2
Accélération moyenne modifiée (1 + γ )2 / 4 1/ 2
u ( i + 1) (i )
= u + ∆t uɺ (i ) (i )
+
( ∆t ) (i ) 2
uɺɺ(i )
2 (2.51)
(i )
∆
uɺ (i +1) = uɺ (i ) +
t
2
(
uɺɺ(i ) + uɺɺ(i +1) )
Bien que l’accélération uɺɺ(i + 1) intervienne dans le calcul de la vitesse, il est possible de calculer le
déplacement à l’incrément i , c’est pourquoi il s’agit d’un schéma explicite. En exprimant
l’accélération dans (2.51)a et en substituant dans (2.51)b, les relations suivantes sont obtenues (pour
un pas de temps constant) :
u (i +1) − u (i −1)
uɺ (i ) =
2 ∆t
( i + 1)
(2.52)
u − 2u (i ) + u (i −1)
uɺɺ(i ) =
( ∆t )2
u (i +1) − u (i )
uɺ (i +1/ 2 ) =
∆t (i +1/ 2 )
uɺ (i +1/ 2 ) − uɺ (i −1/ 2 ) (2.53)
uɺɺ(i ) =
1
2
(
∆t (i + 1/ 2 ) + ∆t (i −1/ 2 ) )
Le principal intérêt d’écrire le schéma sous cette forme est de faciliter son implémentation. Les
déplacements et les vitesses se mettent alors sous la forme :
A chaque pas de temps, l’accélération uɺɺ(i ) est calculée à l’aide de l’équation de la dynamique
(2.47). Ce type de schéma est donc plus aisé à implémenter qu’un schéma implicite, car il n’est pas
nécessaire de mettre en place un algorithme faisant converger la solution à chaque pas de temps. En
contrepartie, la stabilité de ce schéma explicite est conditionnée par le pas de temps : celui-ci doit
(i )
rester en dessous d’une valeur critique ∆t crit afin que la stabilité du schéma soit assurée. Ce pas de
temps critique correspond au temps de propagation d’une onde de compression dans un élément. Il
dépend donc de la plus petite dimension d’élément dans la discrétisation en éléments finis, de la
densité du matériau et de sa rigidité. Le matériau pouvant avoir un comportement non-linéaire,
l’incrément de temps critique doit être réévalué à chaque pas de temps. Un coefficient de sécurité α
est généralement mis en place afin de s’assurer que le pas de temps critique n’est pas dépassé :
∆t (i ) = α ∆t crit
(i )
(2.55)
Abaqus/Explicit utilise ce schéma explicite aux différences centrées à pas de temps adaptatif
[ABA08]. Un organigramme de principe de l’implémentation de ce schéma est proposé sur la figure
2-2. Le calcul du pas de temps critique associé à chaque élément est effectué par Abaqus de la façon
suivante :
ρ
(i )
∆t crit = min Le e (2.56)
éléments λ + 2 µe
pour les matériaux faiblement non-linéaires et/ou faiblement anisotropes. Pour les matériaux étudiés
dans ce manuscrit, au comportement fortement anisotrope, un pas de temps critique maximum sera
calculé manuellement et imposé à Abaqus, à l’aide de la formule suivante :
ρ
(i )
∆t crit = min ( Le ) (2.57)
E éléments
où E est une estimation du module d’élasticité maximal du matériau, pour le calcul considéré.
Pour une mèche de renfort tissé par exemple, en notant K elong la rigidité maximale de la mèche dans la
direction des fibres et S0 sa section dans la configuration initiale, le module E sera défini par :
K elong
E= (2.58)
0.8 S0
u ( 0 ) = 0 uɺ ( −1 / 2 ) = uɺ ( 0 ) = uɺ0
INITIALISATION
σ ( 0) = σ 0 n = 0 t = 0
u Calcul et concentration de la
matrice de masse M
u (i +1)
ext , ( i )
Calcul des forces externes f
u (i )
Calcul des forces internes
u (i −1) Calcul des différentes mesures de
F PI , ( i )
Figure 2-2. Schéma de principe de l'implémentation d'un code de calcul explicite en différences
centrées à pas adaptatif
Le coefficient 0.8 permet de prendre en compte une éventuelle diminution de la section transverse
de la mèche au cours du calcul, induisant une modification du module apparent. La figure 2-3 donne
une estimation du pas de temps critique pour une mèche de fibres de verre, en fonction de la taille
caractéristique de l’élément. L’algorithme de contact par pénalité étant susceptible de modifier
légèrement le pas de temps critique (une diminution maximale de 4% est indiquée dans le manuel
Abaqus), un coefficient de sécurité sera appliqué sur cette estimation lors des calculs.
1.20E-07
1.00E-07
8.00E-08
∆t crit (s)
6.00E-08
4.00E-08
2.00E-08
0.00E+00
Ainsi que la figure 2-2 le montre, la loi de comportement (et donc la subroutine VUMAT)
intervient dans le calcul des forces internes, lorsque les contraintes σ PI , (i ) induites par les
(i )
déplacements u sont calculées. Les détails de l’implémentation de la subroutine VUMAT sont
présentés dans la partie suivante.
G2 F g
2
γ
O2 e2 = R ⋅ O 2 λ
O1 g
1
G1 u e1 = R ⋅ O 1
α β
a. b.
Figure 2-4. Cas-test utilisé pour la vérification des orientations dans Abaqus. (a) Etat initial. (b)
Transformation de corps rigide ( u , β ) + cisaillement ( γ ) + élongation ( λ )
Les projections théoriques des tenseurs F et U seront d’abord déterminées puis comparées à celle
{
obtenues dans Abaqus dans le cas particulier α = 20° ; β = 15° ; λ = 1.5; γ = 1.2 . Afin de simplifier
les formules, le cosinus et le sinus de l’angle α seront respectivement notés cα et sα , et ceux de
}
l’angle β seront respectivement notés cβ et sβ .
⌢
F = g ⊗ G i = F ij G i ⊗ G j = F ij O i ⊗ O j = Fˆ ij ei ⊗ e j (2.59)
i
G i = G i , O i = O i , ei = ei (2.60)
(
F ij = G i ⋅ F ⋅ G j = G i ⋅ g ⊗ G k ⋅ G j
k
)
( k
)
= G i ⋅ g ⋅ (G k ⋅ G j ) = δ kj G i ⋅ g
k
(2.61)
= Gi ⋅ g
j
λc 0
β − sβ + γ cβ
g = λ sβ , g = cβ + γ sβ , g = 0 (2.62)
1 G 2 G 3 G
0 0 1
On a donc :
⌢
F ij = O i ⋅ F ⋅ O j = O i ⋅ ( F kl G k ⊗ G l ) ⋅ O j
(2.64)
= F kl (O i ⋅ G k ) (G l ⋅ O j )
On introduit le tenseur ∆ de rotation de la base G à la base O , ainsi que ses composantes dans la
base G :
∆ = ∆ ij G i ⊗ G j et O i = ∆ ⋅ G i (2.65)
(( ∆ ⋅ G ) ⋅G ) (G ⋅ ( ∆ ⋅G ))
⌢
F ij = F kl i k l j
(2.66)
= (∆ )
ik
T
F ∆
kl lj
λ cα2 cβ + ( cβ + γ sβ ) sα2 (γ c − s ) c − λ s s
β β
2
α
2
α β
0
( β β ) α α
+ γ c + ( λ − 1) s c s + ( (1 − λ ) c + γ s ) c
β β α sα
2
λ cα sβ + ( −γ cβ + sβ ) sα (c + γ s ) c + λs c
2 2 2
β β α α β
F = 0
+ ( ( 1 − λ ) cβ + γ sβ ) cα sα − ( γ c + ( λ − 1) s ) c
O
β β α sα
0 0 1 (2.67)
1.843 0.694 0
= 0.182 0.883 0
α = 20° , β =15°
λ = 1.5, γ =1.2
0 0 1
Fˆ ij = ei ⋅ F ⋅ e j
⌢ (2.68)
= F kl ( ei ⋅ O k ) (O l ⋅ e j )
Fˆ ij = ei ⋅ F ⋅ e j
(( )( ( ))
⌢
)
= F kl R ⋅ O i ⋅ O k O l ⋅ R ⋅ O j
(2.70)
) (O ⋅ ( R (O ) ))
⌢
= F kl ( ( R (O
mn m ⊗ On ) ⋅ Oi ) ⋅ O k l op o ⊗ O p ⋅O j
⌢
= R ikT F kl R lj
Finalement,
1.651 0.839 0
F = 0.326 1.074 0 (2.71)
E α = 20° , β = 15°
λ = 1.5, γ = 1.2
0 0 1
Analyse du résultat
La déformation précédente est reproduite à l’identique dans un calcul Abaqus (figure 2-5). Les
valeurs du gradient de la transformation et des dilatations, fournies par Abaqus à la subroutine, sont
stockées à chaque incrément.
Figure 2-5. Implémentation du test des orientations dans Abaqus. La base représentée est la base
de Green-Naghdi présentée précédemment.
Après calcul, les projections des tenseurs représentatifs de la transformation fournis par Abaqus
sont les suivantes :
1.843 0.6942 0
F = 0.1821 0.8825 0
(2.72)
Abaqus
0 0 1
et
1.643 0.1572 0
U = 0.1572 0.9282 0
(2.73)
Abaqus
0 0 1
En comparant les valeurs obtenues avec les valeurs calculées précédemment, on peut conclure sur
la problématique posée : contrairement à ce qu’indique la documentation Abaqus, les composantes du
tenseur gradient de la transformation F et du tenseur des dilatations U sont fournies à la subroutine
VUMAT dans la base de l’orientation initiale de l’élément O 4. Comme expliqué précédemment, cette
base est une base spatiale fixe par rapport au repère global du modèle. Notons toutefois que cette
conclusion n’est valable que lorsqu’une orientation a été attribuée au matériau par le mot clé
*ORIENTATION : dans le cas contraire, la base de projection de F est la base mixte G ⊗ O .
1 T
σ = F ⋅ S ⋅ F
E J E E E
(2.74)
1 T
σ = F ⋅ S ⋅ F
O J O O O
T 2
C = F ⋅ F = U (2.75)
O O O O
4
Ce résultat diffère lorsque le mot clé *ORIENTATION n’est pas utilisé dans la mise en donnée de la
simulation considérée. Ce résultat diffère également dans le cas des subroutines UMAT utilisées pour
implémenter les lois de comportement dans le code de calcul Abaqus/Standard.
Nous avons vu que dans le cadre de l’hyperélasticité, le second tenseur de Piola-Kirchhoff dérive
d’un potentiel énergétique w . Il est donc possible à ce stade de calculer le second tenseur de Piola-
Kirchhoff dans la base de l’orientation initiale :
∂w
S = 2 (2.76)
∂C O
O
T
σ = R ⋅ σ ⋅ R (2.77)
E O O O
1 T T
σ = R ⋅ F ⋅ S ⋅ F ⋅ R
E J O O O O O
1
= R T ⋅ F ⋅ S ⋅ F T ⋅ R (2.78)
J O O O
1
= U ⋅ S ⋅ U T
J O O O
Dans la formule (2.78), toutes les grandeurs du terme de droite sont connues. La structure de la
subroutine VUMAT, résultant des équations précédentes, est décrite sur la figure 2-6.
Initialisation de la subroutine
F i , U i , C i
O O O
Calcul des invariants classiques de
la transformation
I 1i , I 2i , I 3i , I 4i , I 5i
Calcul des contraintes de Piola-
Kirchhoff à l’aide de la loi de
comportement
∂wi
Si = 2 i
O
∂C O
Calcul des contraintes de Cauchy
avec (2.78)
σ i = 1 U i ⋅ S i ⋅ U iT
E J O O O
Fin de la subroutine
Conclusion du chapitre 2
A travers ce chapitre, les bases de la mécanique des milieux continus ont été abordées. La
formulation générale des lois de comportement hyperélastiques isotropes et anisotropes en grandes
déformations a été introduite. Ce formalisme nous permettra, dans les chapitres qui suivent, de définir
de nouvelles lois de comportement adaptées aux renforts tissés.
Le principe de fonctionnement du code de calcul Abaqus/Explicit a également été décrit dans les
grandes lignes. La compréhension du fonctionnement des codes de calcul, industriels ou non, est
primordiale lorsque les calculs réalisés sortent du cadre des matériaux élastiques linéaires en petites
perturbation : elle permet d’avoir une vision globale sur l’algorithme utilisé et sur les différents
problèmes qui peuvent être rencontrés. L’un des objectifs de cette section était donc d’avoir les clés
nécessaires à une utilisation du code de calcul adaptée au problème traité, afin de pouvoir résoudre les
difficultés liées à l’implémentation des lois de comportement hyperélastiques pour les renforts tissés.
Finalement, la méthodologie nécessaire à d’implémentation de lois de comportement dans
Abaqus/Explicit a été décrite. Nous disposons donc du formalisme et des outils nécessaires à la
création et à l’implémentation de lois de comportement hyperélastiques adaptées au comportement
mécanique des renforts tissés de composites.
Sommaire du chapitre
Plan d’isotropie
Direction
des fibres
Figure 3-1. Isotropie transverse de la mèche de renfort tissé. Image obtenue par tomographie X
a. b. c. d.
Figure 3-2. Modes de déformation de la mèche (la flèche désigne la direction des fibres) :
(a) élongation (b) compaction (c) distorsion (d) cisaillement transverse
m = F ⋅M (3.1)
F −T ⋅ N 1 F −T ⋅ N 2
n1 = et n 2 = (3.2)
F −T ⋅ N 1 F −T ⋅ N 2
Afin que ces vecteurs soient orthogonaux et normés en configuration initiale comme en
configuration actuelle, et pour que les plans qu’ils définissent contiennent la direction M
d'anisotropie du matériau, ces vecteurs doivent vérifier les conditions suivantes :
N 1 ⋅N 2 =0, M ⋅N 1 = 0 , M ⋅N 2 =0,
n1 ⋅ n 2 = 0 , m ⋅ n1 = 0 , m ⋅ n 2 = 0 , (3.3)
N 1 = N 2 =1
La condition
N 1 ⋅ C −1 ⋅ N 1 ≤ N 2 ⋅ C −1 ⋅ N 2 (3.4)
sera également requise, qui porte uniquement sur le choix de l’ordre des vecteurs N 1 et N 2 . Il a
été montré dans [CRI01] que les conditions précédentes peuvent être satisfaites en définissant les
vecteurs N 1 et N 2 comme les vecteurs propres de la sous matrice de C −1 obtenue en
{M , N 1 , N 2 }
supprimant la première ligne et la première colonne :
−1 −1
C 22 C 23
{N 1 , N 2 } = VecteursPropres −1 −1
(3.5)
C 23 C 33
{m, n1 , n 2 } = Q ⋅ {M , N 1 , N 2 } (3.6)
Le tenseur gradient de la transformation peut alors être écrit dans un système de coordonnées dans
lequel l'observateur serait fixe par rapport à la base {M , N 1 , N 2 } :
Fɶ = Q T ⋅ F (3.7)
Rappelons que si Fɶ et F ne sont pas identiques, le principe d’objectivité impose que les énergies
de déformation découlant de ces deux tenseurs soient identiques. Par ailleurs, les relations suivantes
−1
montrent que M est un vecteur propre de Fɶ , ainsi que de son inverse Fɶ :
Fɶ ⋅ M = Q T ⋅ F ⋅ M = M ⋅ C ⋅ M Q T ⋅ m = M ⋅ C ⋅ M M (3.8)
−1 −1 −1 1
M = Fɶ ⋅ Fɶ ⋅ M = M ⋅ C ⋅ M Fɶ M ⇒ Fɶ M = M (3.9)
M ⋅C ⋅ M
Les propriétés (3.8) et (3.9) permettent de faire apparaître les relations suivantes :
−1 −1 −1 −1
N 1 ⋅ Fɶ ⋅ M = 0 , N 2 ⋅ Fɶ ⋅ M = 0 , N 1 ⋅ Fɶ ⋅ N 2 = 0 , N 2 ⋅ Fɶ ⋅ N 1 = 0 (3.10)
−1
Ces relations impliquent que la matrice de Fɶ dans la base {M , N 1 , N 2 } est de la forme
suivante :
f −1 f m−11 f m−21
m
Fɶ −1
{M , N 1 , N 2 } = 0 f 11−1 0 (3.11)
0 0 f 22−1
f f m1 fm2
m
Fɶ = 0 f 11 0 (3.12)
{M , N 1 , N 2 }
0 0 f 22
Tout tenseur gradient de la transformation admet donc une représentation sous cette forme. Comme
nous le verrons par la suite, cette représentation permet une interprétation physique simple de chacune
des composantes de Fɶ . Elle permet également, par la forme de la matrice obtenue, de réaliser
simplement des décompositions multiplicatives du gradient de la transformation, et notamment
d’isoler chacun des modes de déformation de la mèche qui sont sollicités.
f
0 0 1 0 0
m
Fɶ = 0 1 0 0 f 11 f 22 0
{M , N 1 , N 2 }
0
0 1 0 0 f 11 f 22
Felong Fcomp
B B
(3.13)
1 0 0 f m1 fm2
1 fm fm
× 0 f 11 / f 22 0 0 1 0
0 0 f 22 / f 11 0 0 1
Fdist Fsh
B B
Ce faisant, les différents modes de déformation de la mèche de renfort tissé sont mis en évidence :
les matrices Felong , Fcomp , Fdist et Fsh obtenues permettent de décrire respectivement l’élongation
dans la direction des fibres, la compaction et la distorsion dans le plan transverse, et le cisaillement
transverse. En notant
f m21 + f m22 f
α elong = f m , α comp = f 11 f 22 , α dist = f 11 / f 22 , α cis = 2
, tan ( γ ) = m 1 (3.14)
fm fm2
α 0 0 1 0 0
elong
Fɶ = 0 1 0 ⋅ 0 α comp 0
{M , N 1 , N 2 }
0 0 1 0 0 α comp
(3.15)
1 0 0 1 α cis cos ( γ ) α cis sin ( γ )
× 0 α dist 0
⋅ 0 1 0
0 0 1/ α dist 0 0 1
Finalement, cinq grandeurs α elong , α comp , α dist , α cis et γ sont nécessaires pour décrire la
transformation de la mèche avec cette approche. Ceci est cohérent avec le fait que cinq invariants sont
nécessaires pour décrire la déformation des matériaux isotropes transverses (§[Link], p. 54). Ces cinq
grandeurs sont caractéristiques des modes de déformation précédents, et constituent des invariants de
la transformation car la décomposition proposée peut être effectuée de façon unique quel que soit le
gradient de la transformation, et quelle que soit la base de projection de celui-ci.
2 1
I 1 = α elong
2
(1 + α cis2 ) + α comp
2
α dist + 2
α dist
2 1 2 1 2
I 2 = α elongα compα dist 1 + α cis + 4 + α comp + α elongα compα cis 2 − α dist cos ( γ )
2 2 2 4 2 2 2
α dist α
dist
I 3 = α elongα comp
2 4
(3.16)
I 4 = α elong
2
I 5 = α elong
4
(1 + α cis2 )
Afin d'obtenir les expressions des grandeurs α elong , α comp , α dist , α cis et γ en fonction des invariants
I 1 , I 2 , I 3 , I 4 et I 5 , il faut inverser le système précédent, ce qui donne :
α elong = I 4
I3
α comp =
I4
2
I 1I 4 − I 5 I I − I5
α dist = + 1 4 − 1 (3.17)
2 I 3I 4 2 I 3I 4
I5
α cis = −1
I 42
I3 2 2 I5 1
I2 + I dist + I 4 I 3 I dist 2+ 4
I4 I 4 I dist
cos2 ( γ ) =
I 2 1
I 4 I 3 52 − 1 I dist − 2
I4 I dist
1
I elong = ln ( I 4 )
2
1 I
I comp = ln 3
4 I4
I 1I 4 − I 5
2 (3.18)
1 I 1I 4 − I 5
I dist = ln + − 1
2 2 I 3I 4 2 I 3I 4
I5
I cis = −1
I 42
Ces quantités décrivent les quatre modes de déformation de la mèche. Une énergie de déformation
sera définie pour décrire le comportement de chacun d’entre eux : ces énergies seront écrites en
fonction de ces invariants. Nous avons vu également que la loi de comportement des matériaux
hyperélastiques fait intervenir la dérivée de ces invariants par rapport au tenseur de Cauchy-Green
droit (équation (2.44), Chapitre 3 p. 67). Ces dérivées donnent les « directions » des contraintes
attribuées à chaque mode de déformation. Les dérivées des invariants (3.18) seront calculées et
étudiées ci-après.
1 2 ∂w T
σ= F ⋅S ⋅ F T = F ⋅ ⋅F
J J ∂C
(3.20)
2 ∂w ∂I elong ∂w ∂I comp ∂w ∂I dist ∂w ∂I cis ∂w ∂ cos ( γ ) T
= F ⋅ + + + + ⋅F
J
∂I elong ∂C ∂I comp ∂C ∂I dist ∂C ∂I cis ∂C ∂ cos ( γ ) ∂C
Dans l’équation (3.20), les termes du type ∂w / ∂I définissent l’amplitude des contraintes associées
à chacun des modes de déformation, tandis que les termes du type ∂I / ∂ C représentent les directions
des contraintes associées à chaque mode de déformation. Ces dérivées, calculées à l’aide des dérivées
des invariants classiques, s’écrivent de la façon suivante :
∂I elong 1
= M
∂C 2I 4
∂I comp
1 1
= C −1 − M
∂C 4 I4
I
∂I dist
(
2 I 4 I − ( I 1 I 4 − I 5 )C −1 + I 1 + 5 M − 2 C ⋅ M + M ⋅ C
I4
) (3.21)
=
∂C 4 ( I 1I 4 − I 5 ) − 4 I 3I 4
2
∂I cis 1 1 2I 5
= (
2 C ⋅ M + M ⋅C − 3 M
∂C 2 I cis I 4 I4
)
Une étude quantitative des dérivées (3.21) est mise en place afin de valider les directions de
contraintes associées à chaque invariant. Pour cela les dérivées des invariants seront calculées en
fonction des composantes du gradient de la transformation (3.12) dans la base {M , N 1 , N 2 } .
1 0 0
2 ∂I elong T 1
F ⋅ ∂C ⋅ F = 0 0 0 (3.22)
J {M , N 1 , N 2 } J
0 0 0
Quel que soit le potentiel associé à cet invariant, les contraintes obtenues seront dirigées
exclusivement par la direction des fibres, ce qui est conforme au comportement souhaité. Le terme
1/ J sera discuté par la suite.
0 0 0
2 ∂I comp T 1
F ⋅ ∂C ⋅ F = 0 1/ 2 0 (3.23)
J {M , N 1 , N 2 } J
0 0 1/ 2
Les contraintes obtenues à l’aide de cet invariant seront donc toujours de type sphérique dans le
plan orthogonal à la direction des fibres.
0 0 0
2 ∂I elong T 1
F ⋅ ∂C ⋅ F = 0 1/ 2 0 (3.24)
J {M , N 1 , N 2 } J
0 0 −1/ 2
Les contraintes découlant de cet invariant sont donc réparties de façon déviatorique entre les deux
directions du plan orthogonal à la direction des fibres. La position de la contrainte négative (ici la
direction 3) est imposée par la condition (3.4).
2 F ⋅ ∂I cis ⋅ F T
=
J ∂C {M , N 1 , N 2 }
f m 1 f 11 f m 2 f 22
0
fm f m21 + f m22 fm f m21 + f m22
(3.25)
1 f m 1 f 11
0 0
J f m f m21 + f m22
f m 2 f 22
0 0
f m f m 1 + f m 2
2 2
Cette expression est plus complexe que les expressions précédentes car elle répartit les contraintes
en fonction de la façon dont chacune des deux directions possibles de cisaillement transverse est
sollicitée.
1
Élongation I elong = ln ( I 4 )
2
1 I3
Compaction I comp = ln
4 I4
I 1I 4 − I 5
2
1 I 1I 4 − I 5
Distorsion I dist = ln + − 1
2 2 I 3I 4 2 I 3 I 4
I5
Cisaillement transverse I cis = −1
I 42
nl
W elong ( )
I elong si I elong ≤ I elong
0
( )
w elong I elong = (3.26)
lin
W elong (I )
elong si I elong > I elong
0
5
Cette non linéarité est répétable : le retour élastique se fait par le même chemin, et des mèches plus rigides
ont été testées avec succès. Ceci montre que cette non linéarité est bien liée au comportement de la mèche et non
à une mauvaise préhension de la mèche dans les mors au cours de l’essai.
180
160
140
120
Effort (N )
100
80 K elong
60
40
20 0
K elong
0
Une réponse en effort linéaire correspond à une énergie de déformation quadratique. La partie de
l’énergie d’élongation w lin
elong
associée à la réponse linéaire de la mèche sera donc écrite comme un
polynôme d’ordre deux. Un polynôme d’ordre trois sera utilisé pour la partie w nlelong . Ces deux parties
seront dont respectivement recherchées sous les formes :
nl
w elong ( )
I elong = a0 + a1 I elong + a2 I elong 2 + a3 I elong 3 (3.27)
lin
w elong ( )
I elong = b0 + b1 I elong + b2 I elong 2 (3.28)
K elong − K elong
0 0
K elong
w nl
elong = 0
I 3
elong + 2
I elong (3.29)
6S0 I elong 2S0
(I )
2
lin
w elong = 2
I elong − 0
I elong I elong + 0
elong (3.30)
2S0 2S0 6S0
Avec la loi de comportement (3.20), le second tenseur de Piola-Kirchhoff est donné par :
0
K elong K elong − K elong
0
I elong + 0
2
I elong si I elong ≤ I elong
0
1 S0 2S0 I elong
S elong = M (3.31)
I4 K elong − K elong
0
K elong
− I 0
elong + I elong si I elong > I 0
elong
2S0 S0
précédemment, la partie linéaire de la densité d’énergie peut être utilisée seule, conduisant à une loi de
comportement à un seul paramètre.
p
K comp I comp si I comp ≤ 0
(
w comp I comp = ) 0 si I comp > 0
(3.32)
Il est supposé que dans la gamme de déformations étudiée (celle de la mise en forme du renfort), la
rigidité en compaction de la mèche n’atteint jamais la rigidité du matériau composant les fibres, c’est à
dire que la mèche ne peut pas être suffisamment compactée pour que tous les vides entre les fibres
soient comblés. Le second tenseur de Piola-Kirchhoff correspondant à la compaction de la mèche
dérive de l’expression (3.32) :
−1 1
( )
p p −1
S comp I comp ≤ 0 = − K comp I comp C − M (3.33)
2 I4
Les paramètres à identifier sont donc K comp et p . Comme il a été mentionné précédemment, le
comportement en compaction ne peut être caractérisé expérimentalement aussi simplement que le
comportement en élongation. Il sera donc fait appel à une méthode inverse, utilisant des essais de
traction équibiaxiale pour identifier les deux paramètres.
Plusieurs raisons ont mené à l’utilisation d’un essai de traction équibiaxiale plutôt que d’un essai
d’écrasement pour identifier le comportement en compaction de la mèche : comme cela a été
mentionné dans le Chapitre 1 (§[Link], p. 33), les essais d’écrasement sont très sensibles au
coefficient de frottement entre les outils et le renfort et à la vitesse de l’essai, et il est difficile d’obtenir
une mesure précise de l’épaisseur du renfort au cours de l’essai. Les essais de traction biaxiale font
intervenir la compaction du renfort de façon prépondérante sans présenter les défauts précédents, c’est
pourquoi ils ont été choisis pour identifier les paramètres de compaction de la mèche.
1
w dist ( I dist ) = K dist I dist 2 (3.34)
2
La constante K dist sera identifiée par méthode inverse, en même temps que les coefficients de la loi
de compaction (3.32). Le second tenseur de Piola-Kirchhoff correspondant à la distorsion de la mèche
dérive de l’expression (3.34) :
I
(
2 I 4 I − ( I 1 I 4 − I 5 )C −1 + I 1 + 5 M − 2 C ⋅ M + M ⋅ C
I4
)
S dist = 2 K dist I dist (3.35)
4 ( I 1 I 4 − I 5 )2 − 4 I 3 I 4
Dans les cas où le dénominateur de l’expression (3.35) s’annule, des simplifications apparaissent et
le tenseur des contraintes reste défini. D’un point de vue théorique, cette expression ne pose donc pas
de problème. D’un point de vue numérique, il faudra toutefois traiter à part le cas de la nullité du
dénominateur lors de l’implémentation de la loi de comportement afin d’éviter d’obtenir des
expressions indéterminées.
Cisaillement
pur de la mèche
Elongation de
la mèche
Comme cela a été décrit au Chapitre 1 (§[Link], p. 25), l’énergie de cisaillement transverse de la
mèche est liée au glissement relatif des fibres dans leur direction privilégiée et à la présence
d’ensimage dans la mèche. Cet ensimage est supposé suffisamment influent pour que le comportement
associé à ce mode de déformation soit élastique. Il est supposé également que ce comportement n’est
pas rigidifiant (raideur constante). L’énergie de déformation est alors écrite sous la forme :
1
w cis ( I cis ) = K cis I cis2 (3.36)
2
K cis 2 (C ⋅ M + M ⋅ C ) − 35 M
1 1 2I
S cis = (3.37)
2 I4 I4
Au cours d’un essai de traction simple sur un renfort tissé, ce mode de déformation intervient
fortement : dans la partie non-linéaire de la courbe de tension uniaxiale, correspondant à la diminution
de l’embuvage du réseau sollicité, seuls les modes d’élongation et de cisaillement travaillent
[KAW73b]. Ensuite, dans la partie linéaire, seule l’élongation travaille (figure 3-4). Une fois le
Sept paramètres doivent être identifiés : trois pour l’élongation, deux pour la compaction, un pour
la distorsion et un pour le cisaillement. Notons que les trois paramètres liés au mode de déformation en
élongation seront très aisés à déterminer. Cette loi de comportement a été implémentée dans le code de
calcul Abaqus/Explicit selon la méthode présentée précédemment (2.5.2, page 60).
La loi de comportement précédente a été implantée dans le code de calcul Abaqus/Explicit avec la
méthodologie définie au Chapitre 2 Des tests élémentaires ont été réalisés afin de vérifier
l’implémentation de la subroutine VUMAT, et de s’assurer de la cohérence du modèle implémenté
avec la loi de comportement hyperélastique formulée pour la mèche de renfort tissé. Ces tests ont été
réalisés sur un élément seul, un hexaèdre à huit nœuds et huit points d’intégration appelé C3D8 dans
Abaqus.
Quatre tests, correspondant à chacun des modes de déformation, sont mis en œuvre pour vérifier
l’implémentation de la subroutine. Au cours de ces tests, deux types de grandeurs sont comparés :
• des énergies de déformation : les énergies de déformation issue du potentiel hyperélastique
(3.38) et calculée par Abaqus à partir des contraintes et des déformations sont comparées.
Cela permet notamment de vérifier que les dérivées des invariants par rapport au tenseur de
Cauchy-Green droit (i.e. les termes ∂I / ∂ C ) ont été calculées correctement.
• des efforts : les efforts caractéristiques obtenus par simulation sont comparés aux efforts
théoriques. Ceci permet de vérifier la bonne implémentation des dérivées du potentiel
énergétique par rapport aux différents invariants (i.e. les termes ∂w / ∂I ).
Ces tests, ainsi que les efforts caractéristiques qui seront comparés, sont introduits sur la figure 3-5.
F F λ 0 0
F = 0 1 0
0 0 1
a.
F
1 0 0
F F F = 0 λ 0
0 0 λ
b.
F
F2
1 0 0
F3 F3
F = 0 λ
0
0 0 1/ λ
c. F2
F
1 γ 0
F = 0 1 0
0 0 1
d. F
La figure 3-6 donne la comparaison des courbes d’énergies et d’efforts obtenues pour chacun de
ces cas tests. Sur chacun des graphes, la différence entre les deux courbes est tracée, dont les valeurs
se lisent sur l’échelle de droite.
Les courbes d’énergie issue du potentiel hyperélastique et calculée par Abaqus à partir des
contraintes et des déformations sont très proches. Cela valide le calcul et l’implémentation des
dérivées des invariants par rapport au tenseur de Cauchy-Green droit.
Pour chacun des cas-tests, les courbes d’efforts obtenues sont également très proches des courbes
théoriques. Cela permet de valider de calcul et l’implémentation des dérivées du potentiel de
déformation par rapport aux différents invariants.
Effort (N )
1
100
5.0E-06
0.5 0 0
-100 0 0.01 0.02 0.03 0.04
0 0.0E+00
-200 -0.002
0 0.01 0.02 0.03 0.04 T emps (s)
a. T emps (s)
Energie (mJ)
E ffort simulé
Effort (N )
D ifférence
-2 E ffort théorique 4.0E-06 0.1 1.0E-06
D ifférence
-3 2.0E-06 0.05 5.0E-07
-4 0.0E+00 0 0.0E+00
0 0.01 0.02 0.03 0.04 0.05 0 0.01 0.02 0.03 0.04
T emps (s)
T emps (s)
0.006 2.0E-06
E nergie Abaqus
0.005
E nergie modèle hyper 1.5E-06
Energie (mJ)
0.004 Différence
0.003 1.0E-06
0.002
5.0E-07
0.001
0 0.0E+00
Différence Différence
Effort (N )
6.0E-08 6.0E-08
-0.005 -0.005
4.0E-08 4.0E-08
-0.01 -0.01
-0.015 2.0E-08 2.0E-08
-0.015
-0.02 0.0E+00 -0.02 0.0E+00
0 0.01 0.02 0.03 0.04 0 0.01 0.02 0.03 0.04
T emps (s) T emps (s)
c.
1.5E-06 3.0E-06
0.3 0
1.0E-06 -0.5 2.0E-06
0.2
-1
5.0E-07 1.0E-06
0.1 -1.5
0 0.0E+00 -2 0.0E+00
Taffetas de verre
Le premier renfort étudié est un taffetas de verre : parmi les différents types d’armures, l’armure de
taffetas présente la géométrie la plus simple, qui peut être décrite de façon très fidèle à la géométrie
réelle par le modèle de Hivet [HIV08]. L’utilisation de ce taffetas pour identifier et valider la loi de
comportement proposée permet ainsi de limiter les imprécisions du calcul liées au modèle
géométrique utilisé. La figure 3-7 présente les caractéristiques géométriques principales de ce taffetas :
il s’agit d’un taffetas de verre quasi équilibré.
Chaîne T rame
Sergé de carbone
Plus onéreux que le taffetas de verre précédent mais présentant des caractéristiques mécaniques
plus élevées ainsi qu’une formabilité plus aisée, le sergé de carbone étudié ici est principalement
utilisé dans des applications aéronautiques. Ses caractéristiques géométriques principales sont décrites
sur la figure 3-8. Comme le taffetas de verre précédent, ce renfort présente une bonne symétrie
chaîne/trame. Il est donc considéré comme étant équilibré dans cette étude.
Chaîne T rame
a. b.
Figure 3-9. Reconstitution du renfort complet par translation d’une maille élémentaire (a) non
déformée et (b) soumise à une transformation homogène
Deux vecteurs, non uniques, sont nécessaires pour décrire la périodicité du renfort. Par souci de
simplicité, ces deux vecteurs seront pris orientés par les directions de chaîne et de trame du renfort, et
de normes égales aux largeurs de la maille élémentaire7. Ils seront notés T c et T t (figure 3-10). Deux
points matériels P1 et P1' du renfort sont appariés s’ils sont images l’un de l’autre par un déplacement
d’un nombre entier de translations T c et T t :
6
Certains essais présentent des effets de bord particulièrement importants (e.g. pour le « picture frame test »
[DUM03a, DUM03b, WIL08, ZHU07]), qui seront discutés lorsque les simulations correspondantes seront
présentées.
7
En configuration de référence, la maille élémentaire est carrée pour un renfort équilibré et rectangulaire
pour un renfort déséquilibré.
Tc
P5 P3 P4
P1 P1' = P1 + 1 × T c
Sc
Tt St
P2 P = P2 + 1 × T c
2
'
Sc ∩ S t
Lorsque l’on isole une maille élémentaire du renfort, seuls les points situés sur les bords ont un
point apparié dans la maille élémentaire (figure 3-10) : les points situés sur chaque côté sont appariés
avec les points du côté opposé, et les nombres nc et nt correspondant sont égaux à 0 , 1 ou −1 .
L’ensemble des points appariés dans la direction chaîne sera noté Sc , et l’ensemble des points appariés
dans la direction trame sera noté St :
{
Sc = P ∈ maille / ∃Q ∈ maille, ∃nc ∈ ℤ *, ∃nt ∈ ℤ,Q = P + ncT c + nt T t } (3.40)
St = {P ∈ maille / ∃Q ∈ maille, ∃n ∈ ℤ, ∃n ∈ ℤ *,Q = P + n T
c t c c +n T }
t t
L’intersection de ces deux ensembles ne contient que les quatre coins de la maille élémentaire.
ϕ(X )
tt
Tt
tc
Tc
ϕm ( X ) ϕl ( X )
Lorsque le renfort est soumis à une transformation homogène, les propriétés de périodicité
précédentes sont conservées, mais les vecteurs T c et T t changent : les nouveaux vecteurs t c et t t
sont les images par cette transformation des vecteurs T c et T t (figure 3-11).
Le champ de déplacement ϕ ( X ) associé à la transformation de la maille élémentaire peut être
décomposé en deux parties : un champ de déplacement homogène macroscopique ϕm ( X ) , et un
champ de déplacement local ϕl ( X ) décrivant les déplacements de chaque point de la maille
élémentaire par rapport aux déplacements macroscopiques :
ϕ ( X ) = ϕm ( X ) + ϕ l ( X ) (3.41)
Dans les simulations qui seront mises en place dans ce chapitre, le champ de déplacement
macroscopique sera imposé à la maille élémentaire (e.g. traction biaxiale, cisaillement) tandis que le
champ de déplacement local fera partie des inconnues du problème.
La périodicité du renfort se traduit par la nécessité que les déplacements locaux des points appariés
soient identiques :
En particulier, les coins de la maille élémentaire étant appariés les uns aux autres, leurs
déplacements locaux doivent être identiques. Un déplacement local identique des coins de la maille
élémentaire correspondant à un mouvement de corps rigide, les mouvements locaux des coins de la
maille élémentaire peuvent être supposés nuls, sans perte de généralité :
ϕ l ( Sc ∩ St ) = 0 (3.43)
ϕ ( P1 ) − ϕ ( P2 ) = ϕm ( P1 ) − ϕm ( P2 ) (3.44)
a. b.
suite de réduire le modèle éléments finis de façon significative, et de simplifier les conditions aux
limites. Ainsi que cela a été discuté dans [BAD08a], l’utilisation d’une maille élémentaire comme
celle de la figure 3-12b, possédant des points matériels sur l’ensemble de sa frontière, permet une
application plus aisée des conditions aux limites et notamment du champ de déplacement
macroscopique ϕm ( X ) . C’est donc ce type de maille élémentaire qui sera utilisée par la suite.
Figure 3-13. Modèles éléments finis utilisés pour les différents cas de chargement sur la maille
élémentaire de taffetas.
Les symétries ce cette maille élémentaire permettent, dans certains cas de chargement, de réduire le
modèle éléments finis utilisé de façon conséquente :
• si les symétries du chargement imposé diffèrent des symétries de la maille élémentaire, ou
si les mèches d’un même réseau interagissent de façon significative au cours du
chargement, un modèle de la maille entière doit être utilisé (figure 3-13b). C’est le cas par
exemple lorsqu’une maille élémentaire de taffetas est sollicitée en flexion ou en
cisaillement dans le plan du renfort.
• si le chargement imposé possède les mêmes symétries que la maille élémentaire et s’il n’y
a pas d’interaction entre les mèches d’un même réseau au cours de ce chargement, le
modèle peut être réduit. Dans le cas d’un taffetas soumis à un chargement de compression
ou de traction uni- ou bi- axiale par exemple, le modèle peut généralement8 être réduit au
huitième de la maille élémentaire (figure 3-13a).
Dans le cas où les symétries du renfort et du chargement permettent de réduire le modèle à une
fraction de la maille élémentaire, les conditions aux limites qui sont appliquées ne sont plus des
conditions de périodicité mais des conditions de symétries. Ces conditions de symétries sont toutefois
compatibles avec les conditions de périodicité précédentes.
8
Certains tissus ont une densité de mèches très importante qui induit, même au repos, un contact entre les
mèches d’un même réseau. Pour ces tissus il ne sera donc pas possible de réduire la taille du modèle éléments
finis.
l’embuvage du renfort (voir §[Link], p. 27) : c’est principalement le modèle géométrique qui
permettra d’obtenir une bonne description de la partie linéaire de la courbe de traction uniaxiale.
Dans ce manuscrit il sera fait usage des modèles développés dans sa thèse par Gilles Hivet [HIV05,
HIV08]. Ces modèles permettent d’obtenir une géométrie proche de la géométrie réelle du renfort en
utilisant des paramètres mesurables sur le renfort (embuvages, largeurs de mèches, épaisseur). Comme
cela a été mentionné au chapitre 1, ils présentent également l’avantage d’éviter les interpénétrations au
niveau des zones de contact entre les mèches. Les deux modèles géométriques utilisés pour décrire les
deux renforts étudiés dans ce manuscrit sont présentés ci-après.
Figure 3-14. Comparaison du modèle géométrique du taffetas de verre issu du modèle de Hivet
avec une coupe obtenue par tomographie X.
1 2 3 4 5
Figure 3-15. Modèle géométrique du sergé de carbone. Comparaison avec des images issues de
tomographie X.
a.
b.
Figure 3-16. Gestion de l'orientation dans Abaqus (a) en intégration complète et (b) en intégration
réduite.
Le choix du nombre d’éléments dans l’épaisseur et dans la largeur de la mèche a peu d’influence
sur les résultats de calcul [GAS00]. Des maillages contenant deux éléments dans l’épaisseur sont donc
utilisés. Par contre l’intégration dans l’élément est importante car l’orientation des fibres est définie à
chaque point d’intégration : dès que l’élément est soumis à des transformations non homogènes il
existe des différences d’orientation des fibres au sein du volume représenté par l’élément (Figure
3-16). Si une intégration réduite est utilisée, ces variations d’orientation au sein de l’élément ne
peuvent être prises en compte, ce qui peut mener à une rigidification de certains modes de
cisaillement. Des éléments cubiques à huit nœuds sont utilisés avec une intégration complète (8 points
d’intégration).
0 100
K elong (N) 11660
40
20
0
nous l’avons vu au paragraphe [Link] (page 27), les modes de déformation les plus sollicités au cours
d’un tel essai sont l’élongation et le cisaillement transverse.
50 1000
Expérience
Ecarts cumulés (N )
40 800
0.3
Effort par mèche (N )
1 600
30
3
20 5 400
10
10 200
0 0
120
Paramètre identifié pour le
Simulation
cisaillement transverse 100
Expérience
Effort par mèche (N )
K cis (MPa) ~3 80
60
40
20
La Figure 3-18 présente les courbes de tension uniaxiale obtenues pour différentes valeurs de
raideur du comportement en cisaillement transverse. La courbe simulée la plus proche de la courbe
expérimentale est celle obtenue pour un coefficient K cis environ égal à 3. Cette identification est
approximative car les courbes issues de l’expérience et de la simulation n’ont pas exactement la même
forme : comme le comportement du renfort en tension uniaxiale est fortement dépendant de la
géométrie de la maille élémentaire (de l’embuvage notamment), cette différence de formes entre les
courbes obtenues peut être attribuée à :
• une description géométrique approximative de la maille élémentaire,
• une mauvaise homogénéité de l’éprouvette utilisée dans l’expérience,
• une description trop pauvre du comportement en cisaillement transverse par la loi de
comportement.
Nous verrons par la suite que l’identification du comportement en cisaillement transverse du sergé
de carbone est plus précise que celle qui est obtenue ici. Il est donc probable que la différence de
forme entre les courbes issues de l’expérience et de la simulation constatée pour le taffetas de verre
soit liée au modèle géométrique ou à l’expérience.
yE(2) − yE(1)
(1) ( Einterp
yEinterp = yE(1) + x − xE(1) ) (3.46)
xE − xE
(2)
L’écart entre les deux courbes est alors mesuré par la méthode des moindres carrés :
R = ∑ ( yEinterp − yS )
2
(3.47)
nS
C’est ce résidu R qui doit être minimisé afin que la courbe simulée s’approche au mieux de la
courbe expérimentale. Pour ce faire, plusieurs algorithmes sont envisageables. Le choix s’est porté sur
l’algorithme de Levenberg-Marquardt, décrit ci-après.
deux algorithmes, tirant partie des avantages de chacun d’entre eux. Dans cet algorithme, le calcul des
nouveaux paramètres se fait à chaque itération à l’aide de la formule suivante :
p i +1 = p i − ( H + λ I ) ⋅ p i
−1
(3.48)
140
120 Expérience
Simulation initiale
Effort par mèche (N)
100 Itération 1
80 Itération 2
Itération 3
60
40
20
10%
Paramètres identifiés pour la
8%
Ecart relatif
compaction et la distorsion
6%
K comp (MPa) 17800
4%
p 7.72
2%
K dist (MPa) 0.0608
0%
Simulation Itération 1 Itération 2 Itération 3
initiale
Les paramètres identifiés par la méthode présentée ici permettent d’obtenir une courbe de tension
biaxiale proche de la courbe expérimentale, et ce en peu d’itérations. L’utilisation de paramètres
initiaux réalistes (i.e. proches de la solution recherché) est nécessaire, afin que l’algorithme ne
converge pas vers un minimum local de la fonction résidu. Pour la mèche de sergé de carbone dont les
paramètres seront identifiés par la suite, les paramètres déterminés ici seront utilisés comme
paramètres de départ.
l'élongation Expérience
200
S0 (mm²) 0.537 Simulation
0
I elong 0.00047
150
Effort (N )
0
K elong (N) 22000
50
De même que pour la mèche de verre, une très bonne corrélation entre l’expérience et la simulation
est obtenue très simplement à l’aide du solveur d’Excel.
45 450
E xpérience
40 400
Kcis=0.1
35 Kcis=1 350
Effort par mèche (N )
Kcis=2.56
Ecarts cumulés (N )
30 Kcis=3 300
25 Kcis=5 250
Kcis=10
20 Kcis=20 200
15 Kcis=50 150
Kcis=2.31
10 100
5 50
0 0
a. b.
0 0.1 0.2 0.3 0.4 0 20 40 60
D éformation (%) Kcis
200
Paramètre identifié pour le
180 E xpérience
cisaillement transverse
160
K cis (MPa) 2.31
Effort par mèche (N )
140 Simulation
120
100
80
60
40
20
c.
0
L’hypothèse précédente, disant que seul le comportement en cisaillement transverse travaille lors
d’un essai de traction uniaxiale, a des limites : par exemple, si une raideur nulle est associée au
comportement en compaction, la section diminuera et le cisaillement transverse ne travaillera pas. Des
valeurs « raisonnables » doivent donc être choisies pour le comportement des autres modes de
déformation afin que cette hypothèse soit valable. La sensibilité des simulations de traction uniaxiale
par rapport au comportement des modes de déformation de compaction et de distorsion a été étudiée.
La figure 3-22 montre que cette sensibilité est très faible, malgré de fortes variations des paramètres de
comportement associés à la compaction et à la distorsion.
200 5
Kcomp=1000
Kcomp=10000 4
150
Ecart (N )
E cart 10000-1000
100 E cart 100000-1000
E cart 1000000-1000 2
50
1
0 0
Kdist=1
Kdist=5
3
Ecart (N )
E cart 0.5-0.1
100 E cart 1-0.1
E cart 5-0.1 2
50
1
0 0
250
Expérience
200 Simulation initiale
Effort par mèche (N)
Itération 1
150 Itération 2
Itération 3
100
50
15%
Paramètres identifiés pour la compaction et la
Ecart relatif
distorsion 10%
K comp (MPa) 17812.9
5%
p 6.82261
200
180 mèche
160 Expérience
140 Simulation
k=0.5
Effort par mèche (N)
120 libre
k=1
100
80
k=2
60
40
20
0
uniaxiale sur la mèche et sur le renfort ainsi que la courbe de tension équibiaxiale ont été utilisées pour
l’identification de la loi de comportement. Les courbes obtenues pour k = 0.5 et k = 2 n’ont pas été
utilisées pour l’identification.
150
k=2
100
50
La courbe obtenue par simulation pour k = 2 est très proche de la courbe expérimentale. Par
contre, pour k = 0.5 , un léger décalage des déformations peut être observé entre la courbe
expérimentale et la courbe simulée. De façon générale dans les essais de tension biaxiale, le point de
départ (i.e. le zéro des déformations) est difficile à obtenir précisément, le tissu pouvant être
légèrement déformé lors de sa manutention. Ceci pourrait expliquer ce décalage constaté entre les
courbes expérimentale et simulée.
Figure 3-26. Déformée de cisaillement de la maille élémentaire de taffetas de verre obtenue par la
simulation.
Les couples de cisaillement, calculés à partir d’un essai de picture frame et à partir de la simulation
de cisaillement de la maille élémentaire, sont présentés sur la figure 3-27.
0.6
Couple surfacique de cisaillement
0.4 Simulation
([Link]/mm²)
0.3
0.2
0.1
0 10 20 30 40 50 60
Angle de cisaillement (deg)
Figure 3-27. Couples de cisaillement obtenus expérimentalement et par simulation pour le taffetas
de verre.
1
Picture frame test [DUM03]
Couple surfacique de cisaillement
Simulation e=0.003
Simulation e=0.0015
0.4
0.2
0 10 20 30 40 50 60
Angle de cisaillement (deg)
L’influence de l’introduction d’une tension biaxiale sur le cisaillement de la maille élémentaire est
très nette sur la figure 3-28. La courbe la plus proche de la courbe expérimentale est celle obtenue
pour une tension e = 0.0015 . On peut supposer que pour l’essai utilisé ici une telle déformation ait été
imposée au renfort, mais même avec cette valeur de déformation initiale le début de la courbe reste
mal décrit. La présence de l’ensimage semble donc avoir une influence déterminante pour cette partie
de l’essai.
Figure 3-29. Déformée de cisaillement pur de la maille élémentaire de sergé de carbone obtenue
par la simulation.
0.6
C ouple surfacique de cisaillement
0.4 Simulation
(N .mm/ mm²)
0.3
0.2
0.1
0 10 20 30 40 50 60
Angle de cisaillement (deg)
identifiées à l’aide d’essais expérimentaux. Finalement la loi de comportement sera utilisée pour
simuler un cas de flexion à 0°/90°, un cas de flexion à 45° et un cas d’emboutissage hémisphérique,
faisant intervenir simultanément les différents mode de déformation. L’étude de ces cas de chargement
révèlera un comportement inattendu des interlocks, dépendant de la courbure, et une première
approche sera proposée pour modéliser ce comportement.
Sommaire du chapitre
a. b. c.
d. e. f.
Figure 4-1. Modes de déformation des renforts interlocks. (a) compression transverse,
(d) cisaillement dans le plan, (b,e) élongations, (c,f) cisaillements transverses.
1
I elong = ln ( I 41 )
(4.1)
2
I elong = ln ( I 42 )
∂I elong
1
1
= M
∂C 2 I 41 1
(4.2)
∂I elong
2
1
= M
∂C 2 I 42 2
1 0 0
2 ∂I
1
1
F ⋅ ⋅FT = 0 0 0
elong
J ∂C {M 1 , M 2 M 3 } J
0 0 0
(4.3)
0 0 0
2 ∂I elong T
2
1
F ⋅ ⋅F = 0 1 0
J ∂C {M 1 , M 2 M 3 } J
0 0 0
Avec les invariants (4.1), les élongations dans le sens chaîne et dans le sens trame donnent bien lieu
à des contraintes orientées exclusivement dans la direction chaîne et dans la direction trame.
G3 g = G3 = 1
3
g
3
Compaction
transverse
λ ≠1
G1 = g
1
Figure 4-2. L'écrasement d’un renfort ne peut être associée à la norme du vecteur covariant
transverse
Dans le modèle hyperélastique de la mèche, la compaction est définie par le changement volume
total divisé par l’allongement dans la direction des fibres, ce qui donne un changement de surface du
plan d’isotropie de la mèche. Par extension, la compaction du renfort interlock sera définie par le
changement de volume total ( J = I 3 ) divisé par les déformations en élongation dans les directions
chaîne et trame (resp. I 41 et I 42 ):
1 I3
I comp = ln (4.4)
2 I 41 I 42
Notons que le changement de volume associé à un cisaillement pur du renfort dans le plan sera pris
en compte par l’invariant (4.4), puisque dans le cas d’un cisaillement plan il n’y a pas d’allongement
des directions chaîne et trame. Cet aspect sera discuté plus en détail dans la partie concernant
∂I comp ∂ 1 1 1
= ln ( I 3 ) − ln ( I 41 ) − ln ( I 42 )
∂C ∂C 2 2 2
(4.5)
1 1 1
= C −1 − G1 − G2
2 I 41 I 42
0 0 0
2 ∂I comp T 1
F ⋅ ∂C ⋅ F = 0 0 0 (4.6)
J { M 1 , M 2 M 3 } J
0 0 1
Cette représentation est compatible avec les directions de contrainte recherchées : les contraintes de
compaction obtenues en dérivant l’invariant I comp sont orthogonales au plan du renfort.
γ = ( M 1 , M 2 ) − (m 1 , m 2 ) (4.7)
M2 m2
γ
2
m1
γ
M1
2
Les deux réseaux étant orthogonaux en configuration non déformée, le sinus de cette variation
d’angle est :
π
sin ( γ ) = sin − ( m 1 , m 2 ) = cos ( m 1 , m 2 ) (4.8)
2
Le cosinus des vecteurs covariants m 1 et m 2 s’écrit simplement à partir de leur produit scalaire, ce
qui donne, par substitution dans (4.8) :
m1 ⋅ m 2
sin ( γ ) = (4.9)
m1 m 2
Ainsi que nous l’avons vu précédemment, les normes de ces vecteurs sont données par les
invariants I 41 et I 42 : m 1 m 2 = I 41 I 42 . Par ailleurs, leur produit scalaire peut se réécrire :
m1 ⋅ m 2 = F ⋅ M 1 ⋅ F ⋅ M 2
= Fij M 1 j Fik M 2 k (4.10)
= M 1 ⋅C ⋅ M 2 = M 2 ⋅C ⋅ M 1
I 412 I 421
L’invariant utilisé pour décrire le cisaillement dans le plan sera donc le suivant :
I 421
I cp = = sin ( γ ) (4.11)
I 41 I 42
∂I cp 1 ∂I 421 I ∂I 41 I ∂I 42
= − ct − ct (4.12)
∂C I 41 I 42 ∂C 2 I 41 ∂C 2 I 42 ∂C
∂ ∂A
∂B
(
x⋅A⋅y =x⋅
∂B
⋅y ) (4.13)
∂I 421 ∂
∂C
=
∂C
(
M 2 ⋅C ⋅ M 1 )
∂C
=M2⋅ ⋅M 1 (4.14)
∂C
= M 2 ⋅ I + I T ⋅ M 1
1
2
9
La dérivée d’un tenseur d’ordre 2 symétrique par rapport à lui-même n’est pas intuitive : on pourrait dire
grossièrement que cela vient du fait qu’il existe un lien entre les composantes d’un côté et de l’autre de la
( ) ( )
diagonale. La formule ∂C / ∂C = 1/ 2 ⋅ I + I = 1/ 2 ⋅ δ ij δ kl + δ ikδ jl ei ⊗ e j ⊗ ek ⊗ el peut être trouvée
T
Afin de simplifier l’expression (4.14), les différents tenseurs sont projetés dans une base
orthonormée {ei } :
∂I 421
= ( M 2m em ) ⋅ (δ ij δ kl + δ ikδ jl ) ei ⊗ e j ⊗ ek ⊗ el ⋅ ( M 1n en )
1
∂C 2
=
1
2 ( )
( M 2m em ) ⋅ δ ij δ kl ei ⊗ e j ⊗ ek ⊗ el ⋅ ( M 1n en )
1
( )
+ ( M 2m em ) ⋅ δ ikδ jl ei ⊗ e j ⊗ ek ⊗ el ⋅ ( M 1n en )
2
(4.15)
1 1
= ( M 2i ei ) ⊗ ( M 1k ek ) + ( M 1i ei ) ⊗ ( M 2 k ek )
2 2
1
= (M 2 ⊗ M 1 + M 1 ⊗ M 2 )
2
∂I cp 1 ∂I 421 I cp ∂I 41 I cp ∂I 42
= − −
∂C I 41 I 42 ∂C 2 I 41 ∂C 2 I 42 ∂C
(4.16)
1 I I
= ( M 1 ⊗ M 2 + M 2 ⊗ M 1 ) − cp M 1 − cp M 2
2 I 41 I 42 2 I 41 2 I 42
I 4i 3
I cti = = sin (α i 3 ) ( i = 1, 2 ) (4.17)
I 4 i I 43
Et leurs dérivées :
∂I cti 1 I I
= (G i ⊗ G 3 + G 3 ⊗ G i ) − ct G i − ct G 3 (i = 1, 2 ) (4.18)
∂C 2 I 4 i I 43 2 I 4i 2 I 43
1 I3
Compaction I comp = ln
2 I 41 I 42
I 421
Cisaillement dans le plan I cp =
I 41 I 42
I 413
Cisaillement transverse chaîne I ct 1 =
I 41 I 43
I 423
Cisaillement transverse trame I ct 2 =
I 42 I 43
Tableau 5. Invariants associés aux modes de déformation des renforts interlocks épais.
Ces invariants étant définis, des densités d’énergies de déformation peuvent maintenant leur être
associées afin de décrire le comportement mécanique de chacun des modes de déformation du renfort.
Direction Chaîne
Direction Trame
A l’aide de ce renfort, différents essais de caractérisation expérimentale des renforts interlock épais
ont été développés au laboratoire dans le cadre de la thèse de J.G. Orliac. Ceux qui ont été utilisés
seront présentés au fur et à mesure de l’identification du comportement mécanique des différents
modes de déformation. A l’aide de ces essais, des densités d’énergie de déformation sont attribuées à
chaque mode de déformation (i.e. à chaque invariant). La définition et l’identification des paramètres
de chacune de ces densités d’énergies seront présentées en même temps, bien qu’en pratique la loi de
comportement a été implémentée entièrement avant l’identification.
1 0
( )
2
i
K elong I elong i
si I elong ≤ I elong
0
2
w i
= (4.19)
( )
elong
1 2 1 0 0 i
i
K elong I elong − I elong
0
+ K elong I elong I elong si I i
elong >I 0
elong
2 2
4500
4000
Expérience
3500
M odèle
3000
Effort (N )
2500
2000
1500
1000
500
∂w elong
i
0
2 K elong i
I elong si I elong ≤ I elong
0
1
S ielong = 2 = G (4.20)
∂C (
I 4i i 2 K elong I elong − I elong
0
+ K elong
0
) 0
I elong si I elong > I elong
0
0
I elong = 0.0145
0
K elong = 37.85 MPa (4.21)
K elong = 816.33 MPa
La courbe simulée issue de cette identification est présentée sur la figure 4-5. Notons qu’il serait
aisé d’obtenir une courbe d’effort plus fidèle à l’expérience en introduisant une fonction cubique dans
la première partie de la courbe (lorsque I elong ≤ I elong
0
) : c’est ce qui a été fait pour modéliser le
comportement de la mèche en tension dans la loi de comportement définie au chapitre 3. La rigidité de
ce mode de déformation étant très importante par rapport aux rigidités des autres modes de
déformation, ce mode de déformation est très peu sollicité lors de la mise en forme. Cette complexité
n’a donc pas été introduite.
−p
I comp I comp
w comp = K comp 1 − 0 −p − 1 (4.22)
I comp 0
I comp
0
Le coefficient I comp correspond à la valeur que prendrait l’invariant I comp si le renfort était compacté
jusqu’à l’obtention d’une proportion de 100% de carbone. Il s’agit donc d’une valeur limite de
0
l’invariant I comp . Ce coefficient I comp peut être déterminé très simplement en pesant l’éprouvette : en
notant mepr la masse de l’éprouvette de volume initial V0 , et en supposant connue la densité du
carbone constituant les fibres ρcarb = 1.78e − 6 kg / mm 3 , ce coefficient s’exprime :
1 mepr =-1.12
0
I comp = ln (4.23)
2 ρcarb V0
− p −1
p K comp I comp 1 1
S comp = 0 1 − 0 − 1 C −1 − G 41 − G 42 (4.24)
I comp
I comp I 41 I 42
0
Pour l’identification des coefficients K comp et p , le solveur d’Excel est utilisé : ces paramètres sont
optimisés afin que la courbe théorique corresponde à la courbe expérimentale (figure 4-6a). Les
paramètres identifiés sont les suivants :
0
K comp = 7.57e-3 MPa
(4.25)
p = 2.85
Les contraintes de Cauchy issues de la simulation de compaction obtenues avec ces paramètres sont
comparées sur la figure 4-6b. Un bon accord est obtenu entre la courbe expérimentale et la courbe
issue de la simulation.
-0.2
0.03
-0.3
Contrainte
0.02 -0.4 mesurée
0.01 -0.5 Contrainte
-0.6 simulée
0
-0.7
a. -0.7 -0.5 -0.3 -0.1 b.
Icomp Icomp
Figure 4-6. Energie de déformation (a) et contraintes de compaction transverse (b) obtenues
expérimentalement et par simulation
L l
α d
On cherche à extraire l’énergie de déformation par unité de volume initial liée au mode de
déformation de cisaillement du renfort interlock épais à partir d’une courbe effort/déplacement
mesurée sur un bias extension test (figure 4-8).
700
600
Figure 4-8. Courbe de « bias extension test » obtenue pour un interlock 4 mm. (a) courbe entière
(b) zoom sur le début de la courbe.
Au cours d’un bias extension test, l’aire des zones cisaillées et demi cisaillées (resp. en bleu foncé
et en bleu clair sur la figure 4-7) diminue de façon significative. Nous ne disposons pas d’une mesure
de l’évolution de l’épaisseur de l’éprouvette au cours d’un bias extension test. Il est donc difficile de
savoir si cette variation d’aire induit une variation de volume de la préforme. Néanmoins, l’épaisseur
de la préforme ne pouvant pas varier au niveau des mors, une telle variation de volume existe
nécessairement. Il n’est donc pas envisageable de déduire directement le comportement en
cisaillement dans le plan à partir d’un bias extension test, dans la mesure où l’énergie de compaction,
prenant en compte cette variation de volume, ne sera pas nulle. Par ailleurs, les modes de déformation
en élongation dans la direction des mèches risquent également d’être sollicités car leur raideur est
moins grande que celle d’un renfort tissé 2D, et donc l’hypothèse de rigidité infinie, nécessaire pour
obtenir la cinématique de la figure 4-7, n’est pas respectée.
Afin de déterminer l’énergie de cisaillement dans le plan, une méthode inverse par un algorithme
d’optimisation est adoptée (comme pour l’identification de la compaction de la mèche dans le plan
transverse au chapitre précédent). Une forme d’énergie est postulée, dont les paramètres sont recalés
afin de faire correspondre les courbes efforts/déplacements issues de l’expérience et de la simulation.
L’énergie de cisaillement dans le plan est définie par morceaux. Un premier morceau correspond à la
variation de raideur du comportement en cisaillement dans le plan que l’on peut observer sur le début
de la courbe de bias extension test (figure 4-8b) : un zoom sur cette courbe montre que la raideur
mesurée au cours d’un bias extension test est d’abord relativement importante avant de diminuer de
façon significative. Ce comportement peut correspondre à la rupture de l’ensimage. Il est important de
prendre en compte ce comportement, notamment pour deux raisons :
• pour obtenir le palier d’effort que l’on peut observer sur les courbes de bias extension test
• pour éviter les flambements du mode de déformation en cisaillement dans le plan qui peuvent
apparaître au cours des simulations : si la raideur initiale est choisie nulle (ce qui est
généralement obtenu avec une loi puissance ou exponentielle), le cisaillement dans le plan
devient un mode de déformation à énergie nulle en début d’essai.
Un polynôme d’ordre 3 sera utilisé pour décrire cette variation de raideur en début d’essai. Le
second morceau de l’énergie de déformation sera décrit à l’aide d’une loi puissance. Finalement, la
forme postulée pour l’énergie de déformation en cisaillement dans le plan est la suivante :
( )
w cp1 I cp si I cp ≤ I cp0
( )
w cp I cp = (4.26)
w cp2 (I ) cp sinon
avec :
( )
w cp1 I cp = K cp12 I cp2 + K cp13 I cp3
(4.27)
( I ) = K (1 − I )
− p22
w cp2 cp cp 21 cp + w cp 20
Afin que la raideur évolue de façon continue, la fonction w cp doit être de classe C 2 , c'est-à-dire
qu’elle doit être continue et dérivable, que sa dérivée doit être continue et dérivable et que sa dérivée
seconde doit être continue. Ceci induit des relations entre les différents coefficients des fonctions w cp1
et w cp2 . Une fois ces conditions prises en compte, seules trois grandeurs suffisent pour décrire le
comportement en cisaillement dans le plan :
• la raideur initiale : K cp0
• la grandeur définissant la limite entre les deux formes d’énergies (4.27) : I cp0
• la raideur à la limite entre les deux formes d’énergies (4.27) : K cp1
Avec ces trois grandeurs, les coefficients des énergies (4.27) s’écrivent :
K 1 (2 − 3 I cp0 ) − K 0 I cp0
p22 =
K I cp0 ( K 0 + K 1 )
K cp12 = 0
2 K0 + K1 0
( )
1+ p22
K1 − K0 K 21 = I cp 1 − I cp0 (4.28)
K cp13 = 2 p22
6 I cp0
( )
w 20 = w cp1 I cp0 − w cp2 I cp0 ( )
De la même façon que pour l’identification des coefficients de compaction et de distorsion de la loi
de comportement de la mèche, un algorithme de Levenberg-Marquardt est utilisé afin de déterminer
ces trois grandeurs. Les paramètres identifiés pour l’interlock 4mm considéré sont les suivants :
1000
Expérience
800 Simulation initiale (5.00E-03,4.00E-02)
Boucle 1 (4.32E-03,3.08E-02)
Effort (N)
600 Boucle 2 (4.20E-03,2.95E-02)
400
200
0
a.
b. b. 0 20 40 60 80
D éplacement (mm)
50%
Ecart relatif
0%
Simulation Boucle 1 Boucle 2
c.
c.
47.1% 5.7% 1.6%
Figure 4-9. Identification du comportement en cisaillement à partir d’un essai de bias extension
test. (a) Déformée. (b) Courbes effort/déplacement (valeur des paramètres à chaque boucle entre
parenthèses). (c) Ecarts relatif entre les courbes simulées et la courbe expérimentale.
100
Ecart à la simulation de référence (%)
10
0.1
0.01
Figure 4-10. Sensibilité à la vitesse de traction des simulations de bias extension test
identique, inversement proportionnel à la vitesse de déformation. La figure 4-10 donne les écarts
obtenus par rapport à une simulation de référence pour différentes vitesses de déformation.
La figure 4-10 montre que la sensibilité de la simulation du bias extension test à la vitesse de
traction est faible jusqu’à des valeurs de l’ordre de 10 mm/s. L’erreur augmente ensuite rapidement
avec la vitesse. Les simulations réalisées pour l’identification du comportement en cisaillement dans le
plan ont donc été réalisées à 3mm/s.
a.
60 2000
Ecart à la simulation de
50 REFERENC E
1500 2
40
référence
Effort (N )
5
30 1000 10
20 20
500
10 33
0 0
0 20 40 0 20 40 60
b. Largeur des éléments (mm) c. D éplacement (mm)
Figure 4-11. Sensibilité au maillage des simulations de bias extension test lorsque les éléments
sont alignés avec les bords de l'éprouvette.
8 600
Ecart à la simulation de
500 REFERENC E
6 5.051
400 7.856
Effort (N )
référence
4 300 10.102
17.667
200 23.57
2
100 35.3553
0 0
0 20 40 0 20 40 60
Largeur des éléments (mm) D éplacement (mm)
Figure 4-12. Sensibilité au maillage des simulations de bias extension test lorsque les éléments
sont alignés avec les directions des mèches
Des fils témoins sont déposés sur le côté de l’éprouvette, qui permettent de vérifier l’homogénéité
du champ de déplacement au sein de l’éprouvette et ainsi de déterminer le champ de validité de l’essai.
La figure 4-14 présente la comparaison entre l’angle de cisaillement théorique et l’angle de
cisaillement calculé à partir du traitement des images d’essai. A partir d’environ 22°, les fils se
détachent et le calcul de leur orientation par traitement d’images n’est plus pertinent. Une bonne
corrélation est obtenue entre l’angle théorique et l’angle mesuré pour la gamme de déformation
considérée. L’angle théorique est donc utilisé par la suite pour écrire l’énergie de déformation mise en
œuvre au cours de l’essai.
30
Angle calculé
25
20 Angle théorique
Angle (deg)
15
10
0 2 4 6 8
D éplacement du mors (mm)
phases : une première durant laquelle la raideur diminue de façon significative, et une seconde au
cours de laquelle la réponse en effort est quasi linéaire.
10
Effort (N )
6
4
Chaîne
2
Trame
0 2 4 6 8
D éplacement du mors (mm)
Les courbes de cisaillement transverse mesurées nous permettent d’envisager une énergie de
déformation définie par morceaux pour chacune des phases de la déformation. La première phase sera
décrite par une loi puissance avec une puissance comprise entre 1 et 2, permettant d’obtenir pour la
réponse en effort une fonction de type racine. Pour la seconde phase une énergie quadratique sera
utilisée, permettant d’obtenir une réponse en effort linéaire. La forme d’énergie proposée est donc la
suivante :
si I cti ≤ I cti0
pcti
K cti I cti
w cti ( I cti ) = avec 0 < pcti < 2 et i = C ,T (4.30)
K 2 I cti + K 1 I cti + K 0
2
Une fois de plus, les conditions de compatibilité entre les deux morceaux de l’énergie (4.30)
permettent d’introduire des relations entre les différents paramètres introduits :
2 2
K 1 = K cti pcti ( 2 − pcti ) ( I cti0 )
pcti − 1
(4.31)
Trois paramètres doivent donc être identifiés avec cette forme d’énergie : K cti , pcti et I cti0 . Ces
paramètres sont identifiés sur la courbe expérimentale à l’aide du solveur d’Excel. Les valeurs
obtenues au terme de l’identification sont les suivantes :
Les courbes identifiées obtenues avec la description précédente de l’énergie de déformation sont
présentées sur la figure 4-16. La forme de l’énergie de déformation utilisée pour décrire le
comportement du renfort en cisaillement transverse permet d’obtenir des simulations proches des
résultats d’essais.
9 10
8 Expérience Expérience
7 8
Simulation Simulation
6
Effort (N )
Effort (N )
6
5
4
4
3
2 2
1
0 0
a. I cis b. I cis
• la flexion à 0°/90°, c'est-à-dire avec les directions de chaîne et de trame alignées avec les
bords de l’éprouvette. Ce type de sollicitation présente un comportement macroscopique
caractéristique des renforts interlocks épais.
• la flexion à 45°/45°, c'est-à-dire avec les directions de chaîne et de trame tournée de 45°
par rapport aux bords de l’éprouvette. Le comportement du renfort en flexion à 45° est très
différent de celui de la flexion à 0°/90°, ce qui est également un spécificité du
comportement des interlocks.
• l’emboutissage hémisphérique : il s’agit d’une mise en forme complète, dans laquelle le
renfort est embouti par un poinçon hémisphérique.
Ces trois cas de chargement permettront de comparer la loi de comportement à l’expérience, et
mettront en lumière quelques caractéristiques du comportement liées à la microstructure du renfort,
qui s’avèreront impossibles à décrire par une approche de matériau de Cauchy.
[Link] Expérience
Au cours de cet essai, une éprouvette de 200x30x15mm est sollicitée en flexion trois points. La
direction chaîne est alignée avec le grand côté de l’éprouvette, et la direction trame est orientée à 90°
par rapport à ce même côté, c’est pourquoi cet essai de flexion est dit à 0°/90°. Le renfort est posé sur
deux points d’appui fixes écartés de 116mm, puis sollicité au niveau d’un troisième point d’appui situé
au centre de l’éprouvette (figure 4-18). Un déplacement de 60mm du point d’appui central est
effectué, sollicitant l’éprouvette en flexion.
Des fils sont placés sur le côté de l’éprouvette, et une caméra enregistre des images du profil de
l’éprouvette : cela permet d’analyser les déformations de la ligne moyenne et des sections transverses
de l’éprouvette par traitement d’images (figure 4-19, figure 4-20). Les images obtenues pour un
déplacement de l’appui central de 0mm, 30mm et 60mm sont présentées sur la figure 4-18. Ces images
montrent que lors d’un essai de flexion sur un interlock, la cinématique des sections transverses est
particulière : à l’inverse des sections transverses d’une poutre d’Euler-Bernoulli qui restent
orthogonales à la ligne moyenne, celles d’une éprouvette d’interlock restent quasiment verticales. Les
angles entre les sections transverses de l’interlock et la ligne moyenne varient donc fortement. Ces
images permettent un second constat : les zones de l’éprouvette situées à l’extérieur des appuis fixes
sont presque dans le prolongement de la partie centrale de l’éprouvette (un léger angle peut être
mesuré). Ce serait également le cas pour une poutre d’Euler-Bernoulli. Ce qui est remarquable
cependant, c’est que ces zones extérieures sont déformées en cisaillement alors qu’en principe elles ne
sont pas sollicitées : il est donc nécessaire qu’une densité de moment ait été propagée dans cette zone
de l’éprouvette depuis la zone centrale, ce qui est incompatible avec les hypothèses de matériau de
Cauchy. Il s’agit là d’un premier indice poussant à croire que le matériau homogénéisé « interlock »
n’a pas un comportement de matériau de Cauchy : l’énergie de déformation, bien que fortement liée au
premier gradient de la transformation F , est également liée au second gradient de la déformation
∇ F , ce qui n’est pas pris en compte dans le modèle hyperélastique mis en place précédemment.
Effort (N )
3
2
0 10 20 30 40 50 60
D éplacement appui supérieur (mm)
Un capteur placé dans le prolongement de l’appui central permet d’enregistrer les efforts mis en
œuvre lors de l’essai. La courbe d’efforts obtenue est présentée sur la Figure 4-17. Cette courbe
montre que la raideur mesurée diminue au cours de l’essai, et qu’en fin d’essai l’effort lui-même
diminue : les normales aux contacts entre l’éprouvette et les points d’appui fixes tendent à devenir
horizontales, diminuant ainsi l’effort de réaction sur le point d’appui mobile.
Au cours de l’essai de flexion, un effort maximal d’environ 4N est mesuré. L’éprouvette testée ici
pèse 52.5g, ce qui représente un effort de 0.5N, soit environ 12% de l’effort maximal mesuré. Il faut
donc veiller à prendre en compte le poids de l’éprouvette lors des simulations de flexion.
ODB: FLEXION_INTERLOCK15.9mm_trame_POIDS_IdentifieAvecCis.odb Abaqus/Explicit Version 6.8−2 Tue Sep 20 09:36:19 GMT+02:00 2011
[Link] Simulation
Pour la simulation, l’éprouvette complète a été modélisée par des éléments finis, et les appuis ont
été modélisés par des surfaces analytiques rigides. La loi de comportement hyperélastique développée
précédemment est utilisée avec les paramètres issus de l’identification. Afin de satisfaire les
considérations précédentes concernant le poids de l’éprouvette, la gravité est modélisée dans
l’ensemble du modèle. Les efforts et les déformations mis en œuvre au cours de l’essai et de la
simulation seront comparés.
Deux critères seront utilisés pour comparer la cinématique issue de la simulation avec la
cinématique obtenue expérimentalement : la forme de la ligne moyenne et l’orientation des sections
transverses. La comparaison des lignes moyennes est présentée sur la figure 4-19. L’utilisation du
modèle hyperélastique conduit à une ligne moyenne semblable à celle des essais dans la zone centrale
de l’éprouvette, c'est-à-dire entre les deux points d’appui inférieurs. Dans cette zone, la ligne moyenne
obtenue est très proche de celle obtenue par l’expérience. Les rayons de courbure de cette ligne
moyenne au niveau des appuis sont plus faibles dans le cas de la simulation que lors de l’expérience.
10
Coordonnées Y (mm)
-10
0mm Simulation
-30 30mm Simulation
60mm Simulation
0mm Expérience
-50 30mm Expérience
60mm Expérience
-70
Coordonnées X (mm)
Figure 4-19. Comparaison entre expérience et simulation des lignes moyennes de l'éprouvette
d'interlock au cours d'un essai de flexion.
Par ailleurs, une différence importante des lignes moyennes peut être constatée dans les zones
situées à l’extérieur des points d’appui fixes : lors de l’expérience, ces zones sont sensiblement dans le
prolongement de la zone centrale de l’éprouvette, alors que lors de la simulation elles sont proches de
l’horizontale. Ces différences entre la ligne moyenne simulée et la ligne moyenne obtenue
expérimentalement peuvent être attribuées au fait que la loi de comportement proposée ne prend pas
en compte la partie du comportement du matériau dépendant de la courbure du matériau, i.e. du
second gradient de la transformation, qui a été mise en évidence précédemment.
1 3 5 7 9
35 35
Angle de la section (deg)
35 35
Angle de la section (deg)
30 Section 9 (Simu)
Section 9 (Exp)
20
10
0 20 40 60
D éplacement du point d'appui (mm)
Dans l’essai comme dans la simulation, les sections dans la zone centrale ne restent pas
perpendiculaires à la ligne moyenne de l’éprouvette. Il s’agit là d’une caractéristique importante du
comportement des renforts interlocks, mais faut noter que ces sections ne restent pas non plus tout à
fait verticales. De même que la ligne moyenne, les orientations moyennes des fils déposés lors de
l’expérience sur le côté du renfort sont calculées à l’aide d’un algorithme de traitement d’images. Une
comparaison de l’évolution des orientations de certaines sections transverses est proposée sur la figure
4-20. Les orientations obtenues par la simulation à l’extérieur des points d’appui et proche du point
d’appui central (sections 1, 3 et 9) sont en bon accord avec l’expérience. Cependant, dans la zone
située entre les points d’appui fixes (sections 5 et 7) les angles obtenus entre la verticale et les sections
transverses sont plus faibles lors de la simulation que dans l’expérience : les sections tournent moins.
Ceci est encore un effet du fait que la courbure du matériau n’est pas prise en compte dans la loi de
comportement.
7 E xpérience
6 Simulation
5
Effort (N )
0 10 20 30 40 50 60
D éplacement appui supérieur (mm)
Figure 4-21. Comparaison des efforts de flexion sur interlock issus de l'expérience et de la
simulation.
Finalement, les efforts mis en œuvre au cours de l’essai de flexion peuvent être comparés. La figure
4-21 présente les courbes d’effort obtenues expérimentalement et par simulation. Bien que le début de
la courbe montre un assez bon accord entre l’expérience et la simulation, la décroissance de l’effort de
flexion n’est pas mise en évidence par la simulation. Ceci provient du fait que la cinématique issue de
la simulation ne permet pas d’obtenir des normales au contact se rapprochant de l’horizontale au
niveau des points d’appui fixes. Il n’est donc pas étonnant que l’effort ne cesse de croître sur la courbe
issue de la simulation.
Figure 4-22. Déformées de flexion à 45° obtenues par expérience et par simulation.
Malgré ces différences entre l’expérience et la simulation de la flexion à 45°, ce modèle permet une
première compréhension de la cinématique de cet essai, fondamentalement différente de la
cinématique de l’essai de flexion à 0°/90°. Lorsque les mèches sont orientées à 45° par rapport aux
bords de l’éprouvette, des élongations importantes de la partie supérieure et de la partie inférieure de
l’éprouvette sont possibles : c’est notamment la faible raideur en cisaillement dans le plan qui permet à
l’éprouvette de se comprimer dans la partie supérieure et de s’allonger dans la partie inférieure. Cette
raideur en cisaillement dans le plan étant plus faible que la raideur en cisaillement transverse,
l’éprouvette se déforme presque exclusivement en cisaillement dans le plan.
130
Matrice
106
Renfort
Poinçon
100 130
Les vues de dessus des éprouvettes embouties expérimentalement et par simulation permettent de
comparer la cinématique des bords du renfort. La figure 4-25 montre que la cinématique obtenue par
simulation est compatible avec la cinématique obtenue expérimentalement. Quelques plissements
apparaissent néanmoins au cours de la simulation, qui peuvent une fois de plus être attribués au fait
que le comportement lié à la courbure n’est pas inclus dans la loi de comportement.
a. b.
Figure 4-25. Vues de dessus de l'éprouvette emboutie par le poinçon hémishérique : (a) expérience
[DEL09] et (b) simulation.
a.
b.
Figure 4-26. Angles de cisaillements plans (a) et taux de compression (b) dans le renfort interlock
au cours d'un emboutissage hémisphérique, calculés par simulation.
Figure 4-27. Ajout de poutres dans les directions de chaîne et de trame pour modéliser l'influence
du second gradient de la transformation sur le comportement du renfort.
Une simulation de flexion, avec les mêmes caractéristiques géométriques que celui présenté dans la
partie 4.4.1 (p. 127), est réalisée en utilisant ce nouveau modèle. Une comparaison globale des
cinématiques obtenues expérimentalement et par simulation est proposée sur la figure 4-28, La
cinématique de la ligne moyenne est comparée aux valeurs mesurées par traitement d’images sur la
figure 4-29.
Figure 4-28. Comparaison des déformées de flexion 0°/90° obtenues par l’expérience et avec la
superposition du modèle hyperélastique avec des poutres prenant en compte la part de la loi de
comportement liée au second gradient de la transformation.
L’ajout de poutres au modèle, première approche pour la prise en compte de l’influence du second
gradient de la transformation dans la loi de comportement du matériau, permet d’obtenir une
cinématique de flexion plus proche de l’expérience que l’utilisation de la seule loi de comportement
hyperélastique : au-delà des points d’appui extérieurs, les lignes moyennes obtenues prolongent le
reste de l’éprouvette, et une déformation de l’éprouvette est propagée au-delà des points d’appui fixes.
10
Coordonnées Y (mm)
-10
0mm Simulation
-30 30mm Simulation
60mm Simulation
0mm Expérience
-50 30mm Expérience
60mm Expérience
-70
Coordonnées X (mm)
Figure 4-29. Comparaison de l'évolution de la ligne moyenne et des orientations des sections
transverses obtenues expérimentalement et pour une simulation de flexion à 0°/90° avec un
modèle enrichi de poutres.
Toutefois, on constate qu’avec cette approche la raideur en flexion des poutres est trop faible en
début d’essai : cette raideur ne permet pas de redresser suffisamment les bords de l’éprouvette lorsque
le point d’appui central est à mi-course. A l’inverse, cette raideur est trop importante en fin d’essai :
dans la simulation le rayon de courbure obtenu au niveau du point d’appui central est plus important
que lors de l’essai. Ces deux considérations permettent d’imaginer que la raideur attribuée au
comportement du second gradient n’est pas constante : la diminution de cette raideur déductible de la
fin de l’essai de flexion peut être liée à des flambements des mèches à l’échelle mésoscopique, qui ne
sont pas pris en compte par le modèle avec les poutres. Un modèle de poutres au comportement non-
linéaire, ou une « vraie » formulation de la dépendance de la loi de comportement par rapport au
second gradient de la transformation, est donc nécessaire.
Figure 4-30. Comparaison des déformées de flexion à 45° obtenues par l’expérience et avec la
superposition du modèle hyperélastique avec des poutres prenant en compte la part de la loi de
comportement liée au second gradient de la transformation.
L’essai de flexion à 45° a également été simulé avec cette méthode de superposition de poutres
dans les directions de chaîne et de trame avec le modèle de milieu continu proposé. La déformée
obtenue est présentée sur la figure 4-30. Le comportement observé dans la simulation est plus proche
du comportement observé expérimentalement que lorsque le modèle continu est utilisé seul. La
différence de comportement entre la flexion à 0°/90° et la flexion à 45°, caractéristique du
comportement des renforts interlocks, est bien retrouvé par la simulation.
de base des matériaux orthotropes, permettant de décrire chacun des modes de déformation du renfort
interlock : élongation dans la direction chaîne, élongation dans la direction trame, compaction
transverse, cisaillement dans le plan, cisaillement transverse chaîne et cisaillement transverse trame.
Des données expérimentales, obtenues dans le cadre d’un autre travail de thèse au LaMCoS, ont
permis de caractériser chacun de ces modes de déformation, et d’identifier des densités d’énergie de
déformation bien adaptées à leur comportement. La loi de comportement obtenue a été implémentée
dans une subroutine utilisateur VUMAT, permettant ainsi de simuler deux cas de déformation du
renfort : un essai de flexion à 0°/90° et un essai d’emboutissage hémisphérique. L’analyse des
expériences et des simulations de flexion à 0°/90° nous a permis de mettre en évidence un type de
comportement mécanique spécifique aux renforts tissés interlocks : le comportement mécanique de ces
matériaux ne dépend pas que du premier gradient de la transformation, mais également de la courbure
macroscopique du matériau (i.e. du second gradient de la transformation). Les renforts interlocks ne
sont donc pas des matériaux satisfaisant les hypothèses de matériau de Cauchy.
Une première approche a été introduite afin de prendre en compte cette dépendance de la loi de
comportement à la courbure : deux réseaux de poutres ont été ajoutés dans les directions de chaîne et
de trame. Les résultats obtenus à l’aide de cette première approche ont permis de confirmer que la loi
de comportement homogénéisée du matériau interlock est dépendante du second gradient de la
transformation. Toutefois l’approche introduite est trop grossière pour décrire la complexité de cette
dépendance : des phénomènes de flambage de fibres et de mèches apparaissant aux échelles
microscopique et mésoscopique sont certainement responsables d’une non linéarité de ce
comportement dépendant de la courbure, qui ne peut être modélisée par des poutres linéaires.
L’utilisation de poutres non-linéaires ou d’une « vraie » formulation de milieu du second gradient
permettrait d’améliorer cette loi de comportement et d’obtenir une très bonne description du matériau
interlock homogénéisé.
La maîtrise de l’élaboration des pièces en composites à renforts tissés passe par une connaissance
précise de leur comportement à l’échelle mésoscopique et à l’échelle macroscopique. A l’heure
actuelle, le développement de modèles numériques permet d’analyser et de comprendre le
comportement des renforts tissés à ces deux échelles. A terme, les objectifs sont multiples : réduire le
nombre de tests expérimentaux nécessaire à la caractérisation d’un renfort, concevoir numériquement
des renforts adaptés à une application donnée, étudier la formabilité et la perméabilité d’un renfort,
prévoir l’apparition de défauts et, plus généralement, réduire la durée et le coût du cycle de
conception d’une pièce en composite à renfort tissé.
Dans ce contexte, la formulation de lois permettant de décrire le comportement mécanique des
renforts tissés à l’échelle mésoscopique et à l’échelle macroscopique est incontournable. Ces lois
doivent être anisotropes, non linéaires et formulées en grandes transformations. A l’échelle
mésoscopique, les lois de comportement proposées dans la littérature sont hypoélastiques : la difficulté
principale est alors la définition d’une base dans laquelle la dérivée des contraintes est objective. Ces
modèles n’étant pas basés sur des considérations thermodynamiques, ils ne permettent pas une
maîtrise précise de l’énergie de déformation du renfort. À l’échelle macroscopique, les lois de
comportement existantes sont généralement des lois 2D : le renfort est modélisé comme un continuum
2D, en utilisant les éléments finis de coque ou de membrane. Cette hypothèse de modélisation est trop
restrictive pour décrire le comportement de renforts épais tels que les renforts interlocks. L’objectif du
travail de recherche qui a été présenté dans ce manuscrit était d’évaluer les possibilités offertes par les
approches hyperélastiques pour la modélisation du comportement mécanique des renforts tissés de
composites.
Dans ce cadre, une première loi de comportement hyperélastique a été formulée pour décrire le
comportement de la mèche à l’échelle mésoscopique. La déformation de la mèche a été décomposée
en quatre modes de déformation, traduisant les mécanismes de déformation de la mèche. Des
invariants ont été définis décrivant la façon dont ces différents modes sont sollicités. Des fonctions
densité d’énergie de déformation ont été associées à chacun de ces invariants pour modéliser le
comportement mécanique associé à chacun des modes de déformation, définissant ainsi la loi de
comportement. Une méthode a été proposée pour identifier les paramètres de cette loi de
comportement, à partir d’un essai de traction uniaxiale sur la mèche, d’un essai de traction uniaxiale
sur le renfort, et d’un essai de traction biaxiale sur le renfort. Ce dernier essai a par ailleurs été couplé
à un algorithme d’optimisation de Levenberg-Marquardt pour identifier les paramètres du
comportement dans le plan d’isotropie de la mèche. Cette méthode a été appliquée avec succès pour
une mèche de taffetas de verre et pour une mèche de sergé de carbone. Des essais de tension biaxiale
déséquilibrée ont été simulés afin de valider la loi de comportement. Compte tenu des résultats
d’essais expérimentaux, le comportement de la mèche observé au cours de ces simulations s’est avéré
satisfaisant. Des simulations de cisaillement de la maille élémentaire, dans son état libre et soumise à
une sollicitation de tension biaxiale, ont également été mises en place. Le comportement simulé était
conforme au comportement réel du renfort, notamment pour la description de la rigidification en
cisaillement qui apparaît lorsque les mèches de chaîne et de trame entre en contact latéral. Ces cas de
chargement ont donc permis de valider la loi de comportement hyperélastique proposée pour la
description du comportement mécanique de la mèche de renfort tissé.
Une seconde loi de comportement hyperélastique a été proposée, par extension de la première.
Alors que la plupart des lois de comportement macroscopique sont des lois de comportement 2D, une
approche macroscopique 3D a été mise en place, permettant notamment de modéliser le comportement
des renforts épais de type interlock. Cette loi de comportement a été formulée en suivant le même
canevas que pour la formulation de la loi hyperélastique pour la mèche. Les modes de déformation des
renforts interlocks 2.5D ont d’abord été identifiés et décrits à l’aide d’invariants. Des fonctions densité
d’énergie ont alors été associées à ces invariants, dont les paramètres ont été identifiés en utilisant des
résultats d’essais issus d’une autre thèse en cours au LaMCoS. La loi de comportement ainsi obtenue a
été utilisée pour simuler différents cas de déformation macroscopiques : un essai de flexion à 0°/90°
caractéristique des renforts interlock, et un essai d’emboutissage hémisphérique. Un accord
satisfaisant a été obtenu entre les résultats de simulation et les résultats des essais expérimentaux. Il a
toutefois été montré que le comportement mécanique des renforts interlock est sensible à la courbure
de l’éprouvette, et sort donc du cadre des hypothèses de matériau de Cauchy. Cet aspect n’étant pas
pris en compte dans la loi de comportement proposée, les déformées obtenues présentent des
courbures plus importantes que les éprouvettes réelles. Une première approche a été proposée pour
prendre en compte cette influence de la courbure, dans laquelle deux réseaux de poutres ont été
superposés au modèle continu. Cette approche a permis d’obtenir une déformée de flexion plus proche
de la déformée réelle, mais pas tout à fait satisfaisante : la dépendance à la courbure du comportement
du matériau semble être de nature non-linéaire, certainement du fait de flambages de fibres et de
mèches aux échelles inférieures.
Les deux lois de comportement développées au cours de cette thèse ont permis de montrer que les
approches hyperélastiques sont des alternatives intéressantes aux approches hypoélastiques ou
élastiques utilisées dans la littérature : le suivi simultané d’une ou plusieurs directions de fibres est
possible, en grandes transformations, sans avoir à introduire le formalisme des dérivées convectives.
Par ailleurs, la définition de la loi de comportement commençant par une description des différents
modes de déformation permet d’apporter une modélisation plus physique et plus intuitive du
comportement mécanique étudié.
[ABA08] Abaqus 6.8 Analysis User's Manual. Dassault Systemes Simulia Corp., 2008.
[AIM10] Aimène Y., Vidal–Sallé E., Hagège B., Sidoroff F., Boisse P. A hyperelastic approach
for composite reinforcement large deformation analysis. Journal of Composite
Materials, 2010, vol. 44, n°1, 5–26.
[AGO09] Agoras M., Lopez-Pamies O., PonteCastaneda P. A general hyperelastic model for
incompressible fiber-reinforced elastomers. Journal of the Mechanics and Physics of
Solids, 2009, vol. 57, pp. 268–286.
[ALM98] Almeida E.S., Spilker R.L. Finite element formulations for hyperelastic transversely
isotropic biphasic soft tissues. Computer Methods in Applied Mechanics and
Engineering, 1998, vol. 151, pp. 513–538.
[AST02] ASTM, Standard test method for stiffness of fabrics. American Society for Testing and
Materials, 2002, ch. D1388–96(2002).
[BAD08a] Badel P. Analyse mésoscopique du comportement mécanique des renforts tissés de
composites utilisant la tomographie aux rayons X. Thèse de doctorat LaMCoS, Lyon :
INSA de Lyon, 2008, 149 p.
[BAD08b] Badel P., Vidal–Sallé E., Maire E., Boisse P. Simulation and tomography analysis of
textile composite reinforcement deformation at the mesoscopic scale. Composite
Science and Technology, 2008, vol. 68, n°12, pp. 2433–2440.
[BAS09] Basciano C.A., Kleinstreuer C. Invariant-Based Anisotropic Constitutive Models of
the Healthy and Aneurysmal Abdominal Aortic Wall. Journal of Biomechanical
Engineering, ASME, 2009, vol. 131, n°2, 021009.
[BEN07] Ben Boubaker B., Haussy B., Ganghoffer J.F. Discrete models of woven structures.
Macroscopic approach. Composites Part B: Engineering, 2007, vol. 38, pp. 498–505.
[BIC07] Bickerton S., Simacek P., Guglielmi S.E., Advani, S.G. Investigation of draping and
its effects on the mold filling process during manufacturing of a compound curved
composite part. Composites Part A: Applied Science and Manufacturing, 1997, vol. 28,
pp. 801–16.
[BIG05] Bigaud D., Dreano L., Hamelin P. Models of interactions between process,
microstructure and mechanical properties of composite materials–a study of the
interlock layer-to-layer braiding technique. Composite Structures, 2005, vol. 67, n°1,
pp. 99–114.
[BOE87] Boehler J.P. Applications of Tensor Functions in Solid Mechanics. CISM Course no.
292, Springer-Verlag, 1987.
[BOI95] Boisse P., Cherout A., Gelin J.C., Sabhi H. Experimental study and finite element
simulation of a glass fiber fabric shaping process. Polymer Composites, 1995, vol. 16,
n°1, pp. 83–95.
[BOI04] Boisse P. Mise en forme des renforts fibreux de composites. Techniques de l’Ingénieur,
AM 3734, 2004, 10 p.
[BOI06] Boisse P., Zouari B., Daniel J.L. Importance of in-plane shear rigidity in finite element
analyses of woven fabric composite preforming. Composites Part A: Applied Science
and Manufacturing, 2006, vol. 37, pp. 2201–2212.
[BOI11] Boisse P., Hamila N., Vidal-Sallé E., Dumont F. Simulation of wrinkling during textile
composite reinforcement forming. Influence of tensile, in-plane shear and bending
stiffnesses. Composites Science and Technology, 2011, vol. 71, n°5, pp. 683–692.
[BOR03] Borouchaki H., Cherouat A. Drapage géométrique des composites. Comptes Rendus
Mecanique, 2003, vol. 331, 437–442.
[BUE98] Buet–Gautier K. Analyse et simulation du comportement mécanique des renforts
composites tissés. Thèse de doctorat LMMS, Orléans : Université d'Orléans, 1998, 129
p.
[BUE01] Buet–Gautier K., Boisse P. Experimental analysis and modeling of biaxial mechanical
behavior of woven composite reinforcements. Experimental Mechanics, 2001, vol. 41,
pp. 260–269.
[CAI92] Cai Z., Gutowski T. The 3-D deformation behavior of a lubricated fiber bundle.
Journal of composite materials, 1992, vol. 26, 1207–1237.
[CAO08] Cao J., Akkerman R., Boisse P., Chen J., Cheng H. S., de Graaf E. F., Gorczyca J.
L., Harrison P., Hivet G., Launay J., Lee W., Liu L., Lomov S.V., Long A., de
Luycker E., Morestin F., Padvoiskis J., Peng X. Q., Sherwood, J., Stoilova T., Tao
X. M., Verpoest I., Willems A., Wriggers J., Yu T. X., Zhu B. Characterization of
mechanical behavior of woven fabrics: Experimental methods and benchmark results.
Composites Part A: Applied Science and Manufacturing, 2008, vol. 39, n°6, pp. 1037–
1053.
[CIA88] Ciarlet P. Mechanical Elasticity : three dimensional elasticity, Elsevier, 1988.
[COS09] Cosserat E., Cosserat F. Théorie des Corps Déformables. Hermann, Paris, 1909, 226
p., ISBN 2705669205.
[CRI01] Criscione J.C., Douglas A.S., Hunter W.C. Physically based strain invariant set for
materials exhibiting transversely isotropic behavior. Journal of the Mechanics and
Physics of Solids, 2001, vol. 49, pp.871–897.
[DAL08] Dallies E., Berthereau A. Fibres de verre de renforcement. Techniques de l’Ingénieur,
AM 5132, 2008, 19 p.
[DEB06] DeBotton G., Hariton I., Socolsky E. Neo-Hookean fiber-reinforced composites in
finite elasticity. Journal of the Mechanics and Physics of Solids, 2006, vol. 54, pp. 533–
559.
[DEB10] de Bilbao E., Soulat D., Hivet G., Gasser A. Experimental study of bending behaviour
of reinforcements. Experimental Mechanics, 2010, vol. 50, pp. 333–351.
[DEL09] De Luycker E. Simulation et expérimentation en mise en forme de renforts composites
3D interlocks. Thèse de doctorat LaMCoS, Lyon : INSA de Lyon, 2009, 148 p.
[DES05] Desplentere F.S., Lomov S.V., Woerdeman D.L., Verpoest I., Wevers M.,
Bogdanovich A. Micro-CT characterization of variability in 3D textile architecture.
Composite Science and Technology, 2005, vol. 65, pp. 1920–30.
[DIA04] Diani J., Brieu M., Vacherand J.-M., Rezgui A. Directional model for isotropic and
anisotropic hyperelastic rubber-like materials. Mechanics of Materials, 2004, vol. 36,
pp. 313–321.
[DRI10] Dridi S. Essais de caractérisation des structures tissées. Thèse de doctorat LaMCoS,
Lyon : INSA de Lyon, 2010, 125 p.
[DUM03a] Dumont F. Contribution à l'expérimentation et à la modélisation du comportement de
renforts de composites tissés. Thèse de doctorat LMSP/LM2S, Paris : Université de
Paris VI, 2003, 149 p.
142
Cette thèse est accessible à l'adresse : [Link]
© [A. Charmetant], [2011], INSA de Lyon, tous droits réservés
Bibliographie 143
[DUM03b] Dumont F., Hivet G., Rotinat R., Launay J., Boisse P., Vacher P. Mesures de
champs pour des essais de cisaillement sur des renforts tissés. Mécanique & Industries,
2003, vol. 4, pp. 627–635.
[DUP08] Dupupet G. Fibres de carbone. Techniques de l’Ingénieur, AM 5134, 2008, 15 p.
[DUR10] Durville D. Simulation of the mechanical behaviour of woven fabrics at the scale of
fibers. International Journal of Material Forming, 2010, vol. 3, 1241–1251.
[GAS00] Gasser A., Boisse P., Hanklar S. Mechanical behaviour of dry fabric reinforcements.
3D simulations versus biaxial tests. Computational Materials Science, 2000, vol. 17,
pp. 7–20.
[GHO07] Ghoreishi S.R., Cartraud P., Davies P., Messager T. Analytical modeling of synthetic
fiber ropes subjected to axial loads. Part I: A new continuum model for multilayered
fibrous structures. International Journal of Solids and Structure, 2007, vol. 44, pp.
2924–2942.
[GOD09] Godara A., Mezzo L., Luizi F, Warrier A., Lomov S.V., van Vuure A.W.,
Gorbatikh L., Moldenaers P., Verpoest I. Influence of carbon nanotube
reinforcement on the rocessing and the mechanical behaviour of carbon fiber/epoxy
composites. Carbon, 2009, vol. 47, pp. 2914–2923.
[GRU04] Grujicic M., Chittajallu K.M., Walsh S. Effect of shear, compaction and nesting on
permeability of the orthogonal plain-weave fabric preforms. Materials Chemistry and
Physics, 2004, vol. 86, pp.358-369.
[GUO07] Guo Z., Peng X., Moran B. Large deformation response of a hyperelastic fibre
reinforced composite: Theoretical model and numerical validation. Composites Part A:
Applied Science and Manufacturing, 2007, vol. 38, pp. 1842–1851.
[HAG04] Hagège B. Simulation du comportement mécanique des milieux fibreux en grandes
transformations : application aux renforts tricotés. Thèse de doctorat LMSP, Paris :
ENSAM de Paris, 2004, 262 p.
[HAM07] Hamila N. Simulation de la mise en forme des renforts composites mono et multi plis.
Thèse de doctorat LaMCoS, Lyon : INSA de Lyon, 2007, 136 p.
[HAM07b] Hamila N., Boisse P. A meso–macro three node finite element for draping of textile
composite preforms. Applied Composite Materials, 2007, vol. 14, n°4, pp. 235–250.
[HAM09] Hamila N., Boisse P., Sabourin F., Brunet M. A semi-discrete shell finite element for
textile composite reinforcement forming simulation. International Journal for Numerical
Methods in Engineering, 2009, vol. 79, pp. 1443–1466.
[HAN98] Hanklar S. Modélisation mécanique et numérique du comportement des tissus de
fibres. Simulation du comportement mésoscopique de la maille élémentaire. Thèse de
doctorat LM2S, Paris : Université Paris IV, 1998, 156 p.
[HAR04] Harrison P., Clifford M. J., Long A. C. Shear characterisation of viscous woven
textile composites: a comparison between picture frame ans bias extension experiments.
Composites Science and Technology, 2004, vol. 64, n°10–11, pp. 1453–1465.
[HEA01] Heardman E., Lekakou C. In-plane permeability of sheared fabrics. Composites Part
A: Applied Science and Manufacturing, 2001, vol. 32, n°7, pp. 933–940.
[HIR00] Hirokawa S., Tsuruno R. Three-dimensional deformation and stress distribution in an
analytical/computational model of the anterior cruciate ligament. Journal of
Biomechanics, 2000, vol. 33, pp. 1069–1077.
[HIV03] Hivet G. Modélisation mésoscopique pour le comportement bi-axial et la mise en forme
des renforts de composites tissés. Thèse de doctorat LMSP, Orléans : Université
d'Orléans, 2003, 146 p.
[HIV05] Hivet G., Boisse P. Consistent 3D geometrical model of fabric elementary cell.
Application to a meshing preprocessor for 3D finite element analysis. Finite Elements in
Analysis and Design, 2005, vol. 42, n°1, pp. 25–49.
[HIV08] Hivet G., Boisse P. Consistent mesoscopic mechanical behaviour model for woven
composite reinforcements in biaxial tension. Composites Part B: Engineering, 2008, vol.
39, n°2, pp. 345–361.
[HOL00] Holzapfel G.A., Gasser T.C. A New Constitutive Framework for Arterial Wall
Mechanics and a Comparative Study of Material Models. Journal of Elasticity, 2000,
vol. 61, pp. 1–48.
[HUG00] Hughes T.J.R. The Finite Element Method: Linear Static and Dynamic Finite Element
Analysis. Dover Publications, 2000, 672 p., ISBN 0486411818.
[IND08] Indelicato G. Mechanical models for 2D fiber networks and textiles. Thèse de doctorat
LEMTA, Torino : Univertità degli studi di Torino, Institut National Polytechnique de
Lorraine, 2008, 134 p.
[ITS00] Itskov M. On the theory of fourth-order tensors and their applications in computational
mechanics . Computer Methods in Applied Mechanics and Engineering, 2000, vol. 189,
n°2, pp. 419–438.
[JIA99] Jiang W.G., Yao M.S., Walton J.M. A concise finite element model for simple straight
wire rope strand. Internation Journal of Mechanical Sciences, 1999, vol. 41, 143–161.
[KAW73a] Kawabata S., Niwa M., Kawai H. The Finite Deformation Theory of Plain Weave
Fabrics Part I: The biaxial deformation theory. Journal of the Textile Institute., 1973,
vol. 64, n°1, pp. 21–46.
[KAW73b] Kawabata S., Niwa M., Kawai H. The Finite Deformation Theory of Plain Weave
Fabrics Part II: The uniaxial deformation theory. Journal of the Textile Institute., 1973,
vol. 64, n°2, pp. 47–61.
[KAW80] Kawabata S. The standardization and analysis of hand evaluation, The Textile
Machinery Society of Japan, 1980, 97p.
[KHA10] Khan M.A., Mabrouki T., Vidal-Sallé E., Boisse P. Numerical and experimental
analyses of woven composite reinforcement forming using a hypoelastic behaviour.
Application to the double dome benchmark. Journal of Materials Processing
Technology, 2010, vol. 210, pp. 378–388.
[KOI06] Koissin V., Ivanov D.S., Lomov S.V., Verpoest I. Fibre distribution inside yarns of
textile composite: geometrical and FE modelling, in Proceedings of the 8th International
Conference on Textile Composites (TexComp-8), Nottingham, 2006, CD edition.
[KRE09] Krenk S. Non-linear Modeling and Analysis of Solids and Structures. Cambridge
University Press, 2009, 360 p., ISBN: 978-0511601514.
[LAT11] Latil P., Orgéas L., Geindreau C., Dumont P.J.J., Rolland du Roscoat S. Towards
the 3D in situ characterisation of deformation micro-mechanisms within a compressed
bundle of fibres. Composites Science and Technology, 2011, vol. 71, pp. 480–488.
[LAU08] Launay J., Hivet G., Duong A. V., Boisse P. Experimental analysis of the influence of
tensions on in plane shear behaviour of woven composite reinforcements. Composites
Science and Technology, 2008, vol. 68, pp. 506–515.
[LIV64] Livesey R.G., Owen J.D. Cloth stiffness and hysteresis in bending. Journal of the
Textile Institute, 1964, vol. 55, 516–530.
[LOI08] Loix F., Badel P., Orgéas L., Geindreau C., Boisse P. Woven fabric permeability:
From textile deformation to fluid flow mesoscale simulations. Composites Science and
Technology, 2008, vol. 68, pp. 1624–1630.
144
Cette thèse est accessible à l'adresse : [Link]
© [A. Charmetant], [2011], INSA de Lyon, tous droits réservés
Bibliographie 145
[LOM00] Lomov S.V., Gusakov, A. V., Huysmans G., Prodromou A., Verpoest I. Textile
geometry preprocessor for meso–mechanical models of woven composites. Composites
Science and Technology, 2000, vol. 60, n°11, pp. 2083–2095.
[LOM06] Lomov S.V., Verpoest I. Model of shear of woven fabric and parametric description of
shear resistance of glass woven reinforcements. Composites Science and Technology,
2006, vol. 66, pp. 919–933.
[LOM08] Lomov S.V., Boisse P., de Luycker E., Morestin F., Vanclooster K., Vandepitte D.,
Verpoest I., Willems A. Full–field strain measurements in textile deformability studies.
Composites Part A: Applied Science and Manufacturing, 2008, vol. 39, n°8, pp. 1232–
1244.
[LOM11] Lomov S.V., Gorbatikh L., Kotanjac Z., Koissin V., Houlle M., Rochez O.,
Karahan M., Mezzo L., Verpoest I. Compressibility of carbon woven fabrics with
carbon nanotubes/nanofibres grown on the fibres. Composites Science and Technology,
2011, vol. 71, pp. 315–325.
[MAG01] Magno M., Ganghoffer J.–F. Un modèle mésoscopique en grandes perturbations de
matériaux textiles – Application à l’armure toile. Comptes Rendus Académie des
Sciences de Paris, 2001, vol. 329, pp. 149–152.
[MAR63] Marquardt D.W. An algorithm for least squares estimation of nonlinear parameters.
Journal of the Society for Industrial and Applied Mathematics, 1963, vol. 11, pp. 431–
441.
[MIL07] Milani A.S., Nemes J.A., Abeyaratne R.C., Holzapfel G.A. A method for the
approximation of non–uniform fiber misalignment in textile composites using picture
frame test. Composites Part A: Applied Science and Manufacturing, 2007, vol. 38, pp.
1493–1501.
[MIN65] Mindlin R.D. Second gradient of strain and surface–tension in linear elasticity.
International Journal of Solids and Structures, 1965, vol. 1, n°4, pp. 417–438.
[NAW00] Nawrocki A., Labrosse M. A finite element model for simple straight wire rope
strands. Computers and Structures, 2000, vol. 77, 345–359.
[NEW59] Newmark N.M. A method of computation for structural dynamics, Journal of the
Engineering Mechanics Division, ASCE, 1959, vol. 85, n°3, 67–94.
[OGD84] Ogden R.W. Non-Linear Elastic Deformations, Paris : John Wiley, 1984, 544 p
[PAR00] Parnas R.S. Liquid composite molding. Hanser Garner Publications, 2000, 170 p.,
ISBN 1–56990–287–9.
[PEN04] Peng, X., Cao J., Chen J., Xue P., Lussier D., Liu L. Experimental and numerical
analysis on normalization of picture frame tests for composite materials. Composites
Science and Technology, 2004, vol. 64, pp. 11–2.
[PEN05] Peng, X., Cao J. A continuum mechanics-based non-orthogonal constitutive model for
woven composite fabrics. Composites Part A: Applied Science and Manufacturing,
2005, vol. 36, pp. 859–874.
[POT99] Potter K.D. History of the resin transfer moulding for aerospace applications.
Composites: Part A, 1999, vol. 30, pp. 757–765.
[POT06a] Potluri P., Perez Ciurezu D. A., Ramgulam R. B. Measurement of meso–scale shear
deformations for modelling textile composites. Composites Part A: Applied Science and
Manufacturing, 2006, vol. 37, n°2, pp. 303–314.
[POT06b] Potluri P., Parlak I., Ramgulam R., Sagar T.V. Analysis of tow deformations in
textile preforms subjected to forming forces. Composites Science and Technology,
2006, vol. 66, pp.297–305.
[POT08] Potluri P., Sagar T.V. Compaction modelling of textile preforms for composite
structures. Composite Structures, 2008, vol. 86, pp. 177–185.
[QUA94] Quanshui Z., Boehler J. P. Tensor function representations as applied to formulating
constitutive laws for clinotropic materials. Acta Mechanica Sinica, vol. 10, n°4, pp.
336–348.
[RUD97] Rudd C.D., Long, A.C. Liquid moulding technologies: Resin transfer moulding,
structural reaction injection moulding and related processing techniques. Woodhead
Publishing Limited, 1997, 304 p., ISBN 1–85573–242–4.
[RIV48] Rivlin R.S. Large elastic deformations of isotropic materials, Philosophical
Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences,
1948, vol. 241, pp. 379–397.
[SAL05] Salençon J. Mécanique des milieux continus – Tome 1 : Concepts généraux. Ecole
polytechnique (editions), 2005, 370 p., ISBN 2–7302–1245–0.
[SCH92] Schnur D.S., Zabaras N. An inverse method for determining elastic material properties
and a material interface. International Journal of Numerical Methods for Engineering,
1992, vol. 33, pp. 2039–2057.
[SHE07] Sherburn M. Geometric and mechanical modelling of textiles. Thèse de doctorat,
Université de Nottingham, 2007, 271 p.
[SPE07] Spencer A.J.M., Soldatos K.P. Finite deformations of fibre-reinforced elastic solids
with fibre bending stiffness. International Journal of Non-Linear Mechanics, 2007, vol.
42, pp. 355–368.
[SRE06] Sreprateep K., Bohez E.L.J. Computer aided modeling of fiber assemblies. Computer–
Aided Design & Applications, 2006, vol. 3, n°1–4, pp. 367–376.
[TEN07] Ten Thije R.H.W., Akkerman R., Huétink J. Large deformation simulation of
anisotropic material using an updated Lagrangian finite element method. Computer
Methods in Applied Mechanics and Engineering, 2007, vol. 196, pp. 3141-3150.
[VAN91] Van Der Weeën F. Algorithms for draping fabrics on doubly–curved surfaces.
International Journal for Numerical Methods in Engineering, 1991, vol. 31, n°7, pp.
1415–1426.
[VER05] Verpoest I., Lomov S.V. Virtual textile composites software WiseTex: Integration with
micro–mechanical, permeability and structural analysis. Composites Science and
Technology, 2005, vol. 65, n°15–16, pp. 2563–2574.
[WIL08] Willems A., Lomov S.V., Verpoest I., Vandepitte D. Optical strain fields in shear and
tensile testing of textile reinforcements. Composites Science and Technology, 2008,
vol. 68, n°3–4, pp. 807–819.
[YU05] Yu W.R., Zampaloni M., Pourboghrat F., Chung K., Kang T.J. Analysis of flexible
bending behavior of woven preform using non-orthogonal constitutive equation.
Composites Part A: Applied Science and Manufacturing, 2005, vol. 36, pp.839–850.
[YU06] Yu X., Cartwright B., McGuckin D., Ye L., Mai Y.-W. Intra-ply shear locking in
finite element analyses of woven fabric forming processes. Composites Part A: Applied
Science and Manufacturing, 2006, vol. 37, n°5, pp. 790–803.
[ZHA90] Zhang J.M., Rychlewski J. Structural tensors for anisotropic solids, Archives of
Mechanics, 1990, vol. 42, pp. 267–277.
[ZHO04] Zhou G., Sun X., Wang Y. Multi–chain digital element analysis in textile mechanics.
Composites Science and Technology, 2004, vol. 64, pp. 239–244.
[ZHU07a] Zhu B., Yu T. X., Tao X. M. An experimental study of in–plane large shear
deformation of woven fabric composite. Composites Science and Technology, 2007, vol.
67, n°2, pp. 252–261.
146
Cette thèse est accessible à l'adresse : [Link]
© [A. Charmetant], [2011], INSA de Lyon, tous droits réservés
Bibliographie 147
FOLIO ADMINISTRATIF
TITRE :
Approches hyperélastiques pour la modélisation du comportement mécanique de préformes tissées de composites
Ecole doctorale :
Mécanique, énergétique, Génie civil, Acoustique
Spécialité :
Mécanique, Génie mécanique, Génie civil
RESUME :
La simulation des procédés de mise en forme des composites à renforts tissés de type RTM est un enjeu majeur pour les industries de pointe
mettant en œuvre ce type de matériaux. Au cours de ces procédés, la préforme tissée est souvent soumise à des déformations importantes. La
connaissance et la simulation du comportement mécanique de la préforme à l’échelle macroscopique et à l’échelle mésoscopique s’avère
souvent nécessaire pour optimiser la phase de conception de pièces composites formées par de tels procédés.
Une analyse du comportement mésoscopique des préformes tissées de composites est d’abord proposée. Une loi de comportement
hyperélastique isotrope transverse est développée, permettant de décrire le comportement mécanique de chacun des modes de déformation
de la mèche : élongation dans la direction des fibres, compaction et distorsion dans le plan d’isotropie de la mèche, cisaillement le long des
fibres. Une méthodologie est proposée pour identifier les paramètres de cette loi de comportement à l’aide d’essais sur la mèche et sur le
tissu, et une validation par comparaison avec des essais expérimentaux est présentée.
Une analyse du comportement macroscopique des renforts interlocks est ensuite proposée : une loi de comportement hyperélastique
orthotrope est développée et implémentée. Cette loi, extension de la loi de comportement pour la mèche, est également basée sur une
description phénoménologique des modes de déformation de la préforme. Une méthode d’identification des paramètres de cette loi de
comportement est mise en œuvre, utilisant des essais expérimentaux classiques dans le contexte des renforts tissés (tension uniaxiale,
compression, bias extension test, flexion). Cette seconde loi de comportement est validée par comparaison avec des essais d’emboutissage
hémisphérique.
MOTS-CLES :
Renforts tissés, composites, mise en forme, propriétés mécaniques, analyse mésoscopique, analyse macroscopique, lois de comportements
hyperélastiques, grandes transformations, méthode par éléments finis
Directeur de thèse:
BOISSE Philippe
Composition du jury :
BOISSE Philippe, BRIEU Mathias, LADEVEZE Pierre, LOMOV Stepan, MADEO Angela, MAIRE Eric, OTIN Stéphane, SIDOROFF
François, VIDAL-SALLE Emmanuelle