Abichou Hammadi SMZ0110

Simulation de l’emboutissage à froid par une méthode
asymptotique numérique
Hammadi Abichou
To cite this version:

Hammadi Abichou. Simulation de l’emboutissage à froid par une méthode asymptotique numérique.
Mécanique [physics.med-ph]. Université Paul Verlaine - Metz, 2001. Français. �NNT :
2001METZ010S�. �tel-01775454�
HAL Id: tel-01775454

https://hal.univ-lorraine.fr/tel-01775454
Submitted on 24 Apr 2018
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
AVERTISSEMENT
Ce document est le fruit d'un long travail approuvé par le jury de

soutenance et mis à disposition de l'ensemble de la
communauté universitaire élargie.
Il est soumis à la propriété intellectuelle de l'auteur. Ceci

implique une obligation de citation et de référencement lors de
l’utilisation de ce document.
D'autre part, toute contrefaçon, plagiat, reproduction illicite

encourt une poursuite pénale.
Contact : ddoc-theses-contact@univ-lorraine.fr
LIENS
Code de la Propriété Intellectuelle. articles L 122. 4

Code de la Propriété Intellectuelle. articles L 335.2- L 335.10
http://www.cfcopies.com/V2/leg/leg_droi.php
http://www.culture.gouv.fr/culture/infos-pratiques/droits/protection.htm
I
,
INSTITUT SUPERIEURDE GENIE MECANIQUE

ET PRODUCTIQUE
THESE
Présentéeà
U N I V E R S I TDEE M E T Z -
PAR
HammadiABICHOU
Pour I'obtention du gradede :

DOCTEUR de I' UNIVERSITE de METZ
deI'ingénieur
SPÉCIALrTÉ : Sciences
Option: Mécanique
à froid
Simulationde I'emboutissage
par uneMéthodeAsymptotiqueNumérique
Soutenuele 04juillet 2001devantle jury composéde:
J. C. Gelin LMARC & ENSMM Besançon Rapponeur

Y. Q. Guo Universitéde Reims Rapporteur
M. Gueury Univ. H. PoincaréNancY Examinateur
Y. Chevalier ISMCM-CESTIParis Examinateur
H. BenDhia ECAM Examinateur
G. Ferron LPMM - Metz Examinateur
M. Potier-Ferry LPMM - Metz Directeurde Thèse
ISGMP',Ile du saulry' 57M5METZ Cedex01

l,aboratoirede Physiqueet MécaniquedesMatériaux,uMR, 1045,
Mes remerciementsvont en premier lieu ou professeurMichel Potier-
Ferry gui m'ooccueillidqns son équipede rechercheet gui o su m'orienteret
guidermespremiersposdonsle domoinede lo recherche.
J'exprimemes sincàresremerciements
à MessieursJeqn ClaudeJelin et
Ying Qioo GUO pour ovoir ropporter ce trovoil ovec beoucoupd'ottention.por
feurs remarqueset leurssuggestions
ils m'ontpermis d'enrrevoirde nouvelleset
intéressontesperspectivesà ce trovoil.
Je tiens égalementà remercier MessieursM. Gueury,Y. Chevolier, H. Ben

Dhioet G. Ferronde m'qvoirfoit l'honneur
de porticiperoujury.
le remerciemesomieset omisdu L.P.M.M.pour leur oide el leur soutien

donslq réofisotionde ce trovoil et pourtous les bonsmomentspossés ensemble.
,4 ^o fn**n, mzr parenù et a'hnh rno fo*i//n, en Junirin û nn 7ron n...
RE
edl oîHeouEuNn4:fislrAl
N' inv. 5
IooÂo?9
Cote
slb oa[/o
tl gt^.l"-
PLAN DE LA THESE
Introduction générale
CHAPITRE I Etude biblioeranhique 9
Introduction ll
1.1.Description du procédé 11
1.2. Intérêt d'une simulation numérique du procédé T2
1.3. Modélisations mécaniquesdu procédé de I'emboutissage l5
1.3.1.Modèle rigide Plastique 15
1.3.2.Modèle viscoPlastique 16
1.3.3.Modèle suPerPlastique l6
1.3.4.Modèle élastoPlastique t7
1.4. Eléments finis utilisés dans les simulations numériques de I'emboutissage t7
1.4.1.Approche membrane t7
l.4.2.Approches bidimensionnelleset approchessolides l8
1.4.3.Approche coque 18
1.5. Problème de contact et de frottement l9
Introduction 19
1.5.1.Positiondu problème de contact 20
1.5.2.Conditions de contact 2l
f .5.3. Lois de frottements utilisées dans les simulations numériques 22
1.5.3.a.Modèle de Coulomb 22
1.5.3.b.Modèle de Tresca 24
1.5.3.c.Modèle defrottement non local 24
1.5.4.Méthodes de traitement des problèmes de contact et de frottement 25
-1-
1.5.4.a.Méthode de Lemke 26
1.5.4.b.Méthode desMultiplicateurs de Lagrange 27
1.5.4.c.Méthode de Pénalisation 28
1.5.4.d.Méthode de Lagrange augmentée 29
1.5.4.e.Méthode des élémentsd'interfoce 30
1.6. Méthodes de résolution numériques 31
Introduction 31
1.6.1.Méthodesclassiques 32
1.6.2.Méthodes Asymptotiques Numériques 32
1.7.Conclusions 40
et une coque 41
Introduction 42
2.1. Formulation du problème 44
2.1.1.Descriptioncinématiquede la coque 44
2.1.2.Equation d'équilibre 45
2.1.3.Loi de comPortement 46
2.1.4.Etude du problèmede contact 48
2.1.4.a. Conditions de contact 49
2.1.4.b.Régularisation de la loi de contact 50
2.1.4.c. Géométriede contact en 3D 5l
2.2.NléthodeAsymptotique Numérique 53
2.2.1.Méthodesde Perturbation 53
2.2.2.Les approximantsde Padé 56
2.2.3.Discrétisationpar élémentsfinis 56
2.2.4.Procédurede continuation 58
-2
2.3.Conclusions 59
CHAPITIIE 3 Application à la simulation de I'emboutissaeehémisphérioue 60
Introduction 61
3.1. Simulation de I'emboutissagehémisphérique 61
3.1.1.Présentationdu Problème 6l
3.1.2.Traitement des problèmes non-linéaires dans Abaqus 64
3.1.3.Comparaisondesrésultats 67
3.2. Etude de loordre optimal de troncature des séries 69
3.3. Etude de I'influence des paramètres de régularisation 7l
3.3.1.Influencede la régularisationde la loi de comportement 7l
3.3.2.Influencede la régularisationde la loi de contact 72
3.4. Efficacité de I'algorithme 75
3.5. Etude de I'amincissementde la tôle 77
3.6. Etude de la qualité de la solution 79
3.7. Conclusions 81
CHAPITRE 4 Adoucissementde la loi de contact 83
Introduction 83
4.1.Formulationdu problème 83
4.1.1.Conditionsde contactavecinterpénétration 84
4.1.2.Régularisationdesconditionsde contactavecinterpénétration 84
4.1.3.Géométriede contact 86
4.1.4.Techniquede Perturbation 86
hémisphérique
4.2.Applicationà la simulationde I'emboutissage 88
-3-
4.2.1.Gestionde la valeur d'interpénétration 88
4.2.2. Ajustage des paramètres de régularisation 89
4.2.3.Comparaisondes temps de calcul 92
4.2.4.Comparaisonde la qualité de la solution 93
4.3.Conclusions 94
CHAPITRE 5 Généralisation des surfacesdes outils 95
Introduction 96
5.1. Méthodes utiliséespour la représentation de surfacesde contact 96

5.1.1.Maitlage des surfaces en élémentsà géométriesimple 96
5.1.2.Représentationanalytique des surfaces 96
5.1.3.Carreaux de Bézier 97
5.1.3.a. Approximation polynomiale d' ane courbe 97
5.3.2.b. Approxim atio n p olyn omiale d' u n e surfac e 98
5.2. Représentationdes surfacespar les fractions rationnelles 99

5.3. Représentationsdessurfacesdans le cadre de la M.A.N. 101
5.3.1.Régularisation des surfacescomposéesde surfacessimples 101
5.3.2. Traitements des surfacesarbitraires 104
5.4. Conclusions 106
Conclusions généraleset perspectives r07
Référencesbibliographiques 1 09
-4-
I n tro cluctio n générale
Introduction générale
La réalisation par emboutissage, des pièces à parois minces, est couramment

effectuée dans différents secteurs industriels, notamment en construction automobile.
L'emboutissage est un procédé de mise en forme qui consiste à transforrner une tôle fine
plane en une forme non développable, au moyen de pressesmécaniquesou hydrauliques.
Ces dernières années,pour satisfaire les besoins des industriels, en terme de qualité et de
compétitivité, de nombreux travaux ont été consacrés à la simulation numérique de ce
procédéafin d'optimiser et d'assurerla faisabilitédesproduits.
La résolution numérique d'un problème de mise en forme consiste à déterminer les

différentes configurations de la structure au cours de ses déformations,en tenant compte des
conditionsévolutivesde contact ICNRS/I990]. Les méthodesde résolutionles plus connues
actuellementet les plus largement implémentéesdans les codesde calcul par éléments finis
sont les méthodesitérativesbaséessur le processusde prédiction-correction.
Ces méthodesprocèdent en deux étapes; la première consiste à linéariser le problème non
linéaire de départ et donc à prédire une nouvelle solution, tandis que la deuxième étape
consisteà corriger le résidu d'équilibre global de la structurepar itérations successives.
Une bibliographie importante sur ces méthodesest disponible danstous les livres d'éléments
finis traitant les problèmes non linéaires lZienkiewicz et Taylor/9l1 [Crisfield/91]
[Bathe/96].
Les logiciels actuels de simulation numérique utilisent le plus souvent des

algorithmes de résolution de tlpe Newton-Raphson pour résoudre les problèmes non-
linéaires. Les algorithmes implémentésdans ces logiciels de simulation sont d'autant moins
performants lorsque les non-linéarités s'accumulent. C'est notamment le cas de
I'emboutissagepour lequel sont réunies à la fois les non-linéaritésphysiques (plasticité,
éventuellementI'endommagement),des non-linéarités géométriques(grands déplacements,
grandes rotations et grandes déformations) auxquelles il faut ajouter les non-linéarités de
contactet de frottement.Pour ces raisons,la simulation numériquedes procédésde mise en
-5-
forme par éléments finis a pour principal inconvénient le coût de calcul élevé [Bussy et
Vauchez/891.
Ces dernières années, plusieurs travaux ont été consacrésà I'amélioration de la

convergencedes méthodes de résolution et à I'adaptation du pas de calcul [Crisfreld//83]
la résolution des problèmes
[Crisfield//Sl] [Ramm/8l]. L'équipe de Ladevèzetravaille sur
non linéaires en utilisant des méthodes de calcul non incrémentales [Ladevèze/85]
qui à chaque itération , propose une
lLadevèzel91] : il s'agit d'une méthode itérative
approximation des variables en tout point de la structureet sur la totalité de I'incrément.
En parallèle à ces méthodes,nous avons celles de perturbation qui existent depuis
longtemps [Signorinii30]. Ces demières n'avaient pas trouvé une adhésion unanime de la
part des numériciens car elles sont baséessur des développementsanalytiquesde la solution
et leurs applications étaient limitées aux problèmes à géométrie simple. Il a fallu attendre
1968 pour que Thompson et Walker [Thompson et Walker/68] présentent une première
association des méthodes de perturbation et des méthodes d'éléments finis mais les
conclusions quant à I'efficacité de ces méthodes n'étaient pas probantes.En 1990, dans un
article consacréaux problèmes de bifurcations perturbées,Damil et Potier-Ferry ont repris
ces méthodeset, grâce à un choix judicieux de la formulation de base (quadratique),ils ont
pu obtenir des développementsjusqu'à des ordres élevés, ce qui a permis d'améliorer la
solution asymptotiqueen exigeantun temps raisonnablede calcul du secondmembre [Damil
et Potier-Ferryl90l.
La méthode asymptotiquenumérique (M.A.N.) a été fiabilisée quelquesannéesplus
tard par B. Cochelin et ses collaborateurs(1994) en présentantun algorithme permettant de
suivre complètement les branches de solution en appliquant la M.A.N. pas à pas dans une
procédurede continuation. Cette technique de continuation a étêperfectionnéeplus tard par
Elhage Hussein et at lBlhage-husseinet atl}}l qui a proposé une technique de continuation
baséesur les approximantsde Padé.
La M.A.N consiste à chercherla solution d'un problème non linéaire sous forme de
série entière. Elle permet de transformer le problème de départ en une succession de
problèmes linéaires admettant le même opérateur de rigidité tangente.A chaque ordre, les
termesnon linéaires sont reportéssur un secondmembre obtenu de façon optimale.
-6-
La M.A.N. est alors appliquée aux problèmesprésentantdes non-linéaritéstelles que
les coques élastiques en rotations modérées, ou aux indicateurs de bifurcations pour les
problèmes à géométrie non linéaire pour ne citer que quelquesexemples.A partir de cette
avancée considérable, la M.A.N. s'est avérée être une solution alternative aux méthodes
itératives classiques.De plus, pour une classe de problèmes,il est démontré que cette
méthode est rapide (moins de matrices de rigidité à décomposer),fiable (automatiqueavec
des pas adaptatifs) et facile à utiliser. On peut citer à titre d'exemples les problèmes de
coquesen grandesrotations[Ammar/96] [Zahrouni et all99l,la mécaniquedes fluides [Tri
et atl96l lCadou/9l1, la viscoplasticitéen grandesdéformations[Brunelot/99),la plasticité
et le
lBralkatlgTl [Zahrouni et al 198f, le problème de contact [Elhage-Hussein/98]
traitement desproblèmesde dynamique [cohelin et compain/99].
Nous présentonsdans ce travail une modélisation tridimensionnelle du procédé de
I'emboutissageoù sont coupléesplusieurstypes de non-linéarités(géométrique,matérielleet
non-linéaritédue au contact).
Contenu de la thèse
Dans le premier chapitre, on présenteune étude bibliographique sur I'emboutissage

comme procédéde mise en forme, en insistantsur I'aspectnumérique.Ensuite on présente
une bibliographie ou un état de I'art des méthodesasymptotiquesnumériquesqui constituent
le cadre de notre simulation du procédé.
Dans le deuxième chapitre, on présente une modélisation du problème de contact

entre un corps rigide et une coque élastoplastique.Une application de I'algorithme obtenu à
la simulation de I'emboutissagehémisphériquesera présentéedans le chapitre 3. Les
résultats du modèle présentémontreront l'efficacité de I'algorithme en les comparantà ceux
obtenus par le code industriel Abaqus. Une étude détaillée des différents paramètresde
régularisation introduits, pour un choix judicieux de leurs valeurs, seraprésentée.
Une améliorationde la modélisationdes conditions de contact dans le cadre de la

M.A.N. fera I'objet du chapitre4. Dans ce chapitre,on proposeune autrerégularisationde la
loi de contact. On désire par cette nouvelle régularisation,adoucir la loi de contact, qui
-7 -
Introduction eénërale
présente une dureté particulière comparée aux autres types de non-linéarités.

L'adoucissementeffectué consiste à autoriser une interpénétrationentre le corps rigide et la
structure étudiée (tôle déformable). L' algorithme est testé sur le problème étudié dans le
chapitre3.
Une extension de ces algorithmes pour une simulation de I'emboutissage et des

procédés de mise en forme, en généralisant les formes des outils est étudiée au dernier
chapitre. On construit les formes géométriquescomplexes à partir des formes simples. On
procède d'une méthode de régularisationdes formes de surfacesde contact dans certainscas
et d'une méthodede recouvrementdans d'autres cas. Dans les deux cas, une procédurede
recherchede la distanceminimale (point de contact/surfacerigide) est employée.
Les résultats obtenus, montreront la validité des algorithmes et l'efficacité des méthodes
asymptotiquespour le traitementdes problèmesqui combinentdes différentesnon-linéarités,
en les comparantà ceux obtenuspar les méthodesde résolutionclassiques.
-8-
Chapitre I Etude libljpgtsphigug.
CHAPITRE 1
Etude bibliographique
Introduction 11
1.1.Descriptiondu procédé 11
1.2.Intérêt d'une simulation numérique du procédé t2
1.3.Modélisationsmécaniquesdu procédéde I'emboutissage 15

1.3.1.Modèle rigide Plastique 15
1.3.2.Modèle viscoplastique t6
1.3.3.Modèle superPlastique l6
1.3.4.Modèle élastoPlastique t7
1.4. Eléments finis utilisés dans les simulations numériques de I'emboutissage l7

1.4.1.Approche membrane t7
1.4.2.Approchesbidimensionnelleset approchessolides 18
1.4.3.Approche coque 18
1.5. Problème de contact et de frottement t9
Introduction 19
1.5.1.Position du problèmede contact 20

1.5.2.Conditions de contact 2T
1.5.3.Lois de frottements utilisées dans les simulations numériques 22
1.5.3.a.Modèle de Coulomb 22
1.5.3.b.Modèle de Tresca 24
1.5.3.c.Modèle defrottement non local 24
1.5.4.Méthodes de traitement des problèmes de contact et de frottement 25

1.5.4.a.Méthode de Lemke 26
1.5.4.b.Méthode desMultiplicateurs de Lagrange 27
1.5.4.c.Méthode de Pénalisation 28
Chapitre I Etude bibliograPhiaue
1.5.4.d.Méthode de Lagrange augmentée 29

1.5.4.e.Méthode des élémentsd'interface 30
1.6.Méthodes de résolution numériques 31

Introduction 31
1.6.1.Méthodesclassiques 32
1.6.2.Méthodes Asymptotiques Numériques 32
1.7.Conclusions 40
-10-
Introduction
L'emboutissageest un procédéde mise en forme des piècesà froid et quelquesfois à

chaud. Il permet, par déformation plastique des tôles minces et au moyen des outils, d'obtenir
et de formes tridimensionnellesplus ou moins complexes.Les
des piècesnon-développables
tôles à emboutir sont généralementdes tôles laminées d'alliages métalliques ductiles.
L'emboutissage est un procédé qui marque sa présence dans presque toutes les industries
mécaniques.
L'emboutissagea une grande importance industrielle, si I'on se réfère au tonnage des
pièces fabriquées.La faible épaisseurdes tôles utiliséesdans la production,permet un gain
important du poids et une économiede la matière. Les pièces obtenuespar ce procédépeuvent
égalementavoir une granderigidité vues leurs formes et leurs conceptions.C'est un procédé
qui permet donc un bon rapport rigidité-poids, ce qui est très souhaité et recherché dans
l'industrie. C'est grâceà ces avantagesque ce procédétrouve une largeapplication.Il est très
utilisé dans la fabrication du mobilier métallique, de l'électroménager,de l'emballage
métallique, de I'aéronautique et surtout de I'automobile ( il suffit ici de dire que sur un
véhicule I'emboutissageconcerne500 pièces environ ). Un autre avantagetrès important pour
ce procédé : il est très bien adapté au travail en série, un outillage spécial permet d'atteindre
des cadencesde productiond'environ 60 piècespar minute'
1. Description du procédé
La machine à emboutir est constituéede trois élémentsprincipaux : un poinçon, une

matrice et un serre-flan. Le poinçon délimite le contour intérieur et applique la force
nécessaireà la mise en forme de la tôle: force d'emboutissage.Quant à la matrice, elle
délimite le contour extérieur de la pièce. Le serre-flan permet, en pressantla tôle contre la
matrice, de contrôler l'écoulement du métal le long du poinçon et d'éviter le phénomènedu
plissement(figure 1.1)
- 11
Chapitre l Etude bibliographiaue
FouU
Serre
Tôle --+ flanc
Figure (l.l) Description de l'outillage de I'emboutisseuse
La méthodede formagedansce procédédépenddu mode d'action du serre-flan.C'est

un formage par expansion si la partie de la tôle située entre la matrice et le serre-flan est
bloquée. Il est appelé formage par rétreint lorsque cette partie du métal a la liberté de glisser
et de s'écoulerpour alimenterla hauteurde la pièce. Pour contrôlerle glissementde la tôle,
une méthodesouventutilisée dans I'industrie, consisteà freiner latéralementla tôle à I'aide
desjoncs de retenuequi l'obligent à seplier et se déplier.
2. Intérêt d'une simulation numérique du procédé
Les facteursprincipaux dont dépendla réussitede I'opération d'emboutissagesont : le

matériau de |a tôle, les paramètresde I'outillage et les conditions opératoires.Les deux
premiers facteurs constituent des thèmes de recherche majeurs en emboutissage. Les
propriétésmécaniquesdes matériaux qui influent directementsur la réussite de l'opération de
l'emboutissagesont:
- La limite d'élasticité( Re )
- La résistanceà la ruPture ( Rm )
- L'allongementà la ruPture(A)
- Le coefficientd'écrouissage(n)
- Le coefficientd'anisotroPie(r)
qui peuvent être déterminéesau moyen d'un simple essaide traction.
-12-
Chapitre I Etude bibliograqhique
D'autrespropriétéssont souventprisesen comptepour juger I'emboutissabilitéd'une tôle par

exemple:
- la courbe limite de formage ( CLF) ou FDL (Forming Limit Diagram),
- le rapportlimite I'emboutissage,LDR ( Limiting Drawing Ratio ),
- la profondeurmaximalede Godet,
- le profondeur limite de pièce hémisphériqueLDH (Limiting Dome Heith).
Il est très important à ce niveau de tenir compte du fait que ces propriétésne sont pas
constantesau cours de I'opérationde la mise en forme. Ces propriétésvarient en fonction des
cheminsde déformation,la complexitéde la forme despièces,la vitessede déformationet de
l' augmentationde température.
Après avoir déterminé le comportementdu matériau, il faut caractériserI'outillage du procédé
en définissantles propriétéssuivantes:
- la forme et l'état de surfacede I'outillage,
- les profils desoutils,
- le jeu entrele poinçon et la matrice,
- la forme et la position dejonc du retenue,
- la cadencede la production,
- la pressionde serreflan
- la lubrification (dans certain cas)
De nombreux travaux ont été consacrésà l'étude des outils de I'emboutissagepour

améliorerla qualité des produits et réduire les coûts de fabrication.Les outils sont réalisésen
fonte pour des grandesséries ( plus de 1000000de pièces) et en béton de résine pour les
pièces prototypes ( de l'ordre de 100 pièces). Les défauts des pièces emboutiespeuvent
apparaîtrependantou aprèsI'opération de I'emboutissage,ils se produisentsous les formes
suivantes:
- striction /rupture
- plis
- défautdimensionnels( retour élastique,amincissementinacceptable)
- défauts surfaciques(rayures,peau d'orange)
Tous cela montre bien que I'emboutissage,comme procédéde mise en forme, est un
procédé complexe et difficile à maîtriser. Pour cela, des essaisaux laboratoireset aux
-13-
chonitre t Etude bibliosraphique
ateliers sont souventnécessairespour déterminerune combinaisonadéquatede cesparamètres

et obtenir un emboutissageréussi. Or ces essais sont généralement nombreux lents et
coûteux [El Mouatassim et Jameuxlg|] fPicart et al/99]. Ces dernièresannéesdes nombreux
outils de calcul, de simulation et de CAO ont été développés afin de prévoir
I'emboutissabilité des tôles, de réduire les coûts et les délais de production des pièceset
d'améliorer la qualité des produits. Ces outils de calcul et de simulation,qui sont une suite
logique des progrès scientifiquesen calcul des structures,sont devenuspossiblesgrâce à
l'évolution de cette branche de la mécanique et des moyens informatiques à bas prix et aux
développementdes algorithmes bien adaptés.Ces algorithmes sont devenus capablesde
prendre en compte les grandesdéformations,les grandesvitessesde déformation et les lois de
comportementdes matériaux ainsi que le traitement des conditions de contact entre les outils
et la tôle [Afzali et Bouhelier94].
Pour satisfaire les besoins de I'industrie ( et surtout I'industrie automobile) en

simulationnumériquede I'emboutissage,desnombreuxlogiciels ont été développéset mis en
servicepour les industriels.Parmi ces logiciels,on peut citer:
ISOPUNCH(sociétéISOFORMSOLLAC)
SIMEM2(sociétéSIMTECH)
ABAQUS(SociétéABAQUSsoftware)
LS-DYNA-3D(sociétéDYNALYS)
FORGE2(sociétéTRANSVALOR)
OPTRIS(sociétéMATRA DATA VISION)
PAM-STAMP(sociétéESI)
MECALOG)
RADIOSS(société
Tous ces logiciels de simulation sont fondés sur la méthode des élémentsfinis. Dans cette
liste des codeset des outils de simulation du procédé,on peut distinguerdeux catégoriesde
modèles.Des modèlesdits simplifiés et des modèlescomplets.Parmi les modèlessimplifiés,
'approche inverse'. Cette approcheinverse a été
on trouve (par exemple) I'approche dite
développéepar l'équipe de J. L. Batoz lBatoz et all89l. La première version de cette approche
a été présentéeen 1989 en utilisant des élémentsmembranes.Cetteméthodepermet d'éviter
la lenteur de la procédure du calcul incrémentaledans I'analyse de la loi de comportementet
les conditionsde contact.Les points matérielsdans la tôle déforméesont déterminésà partir
-14-
(|hnnitrel - EtudebibliograDhique
de la forme finale de la pièce emboutie(forme CAO maillée). En 1995,cette approchea été

améliorée en tenant compte des phénomènesde flexion qu'on ne peut pas négliger dans la
simulation du procédéfBatozet al/95]. La méthodeinversea été implémentéedans quelques
codesindustrielstels que FAST3D, ICEM-STAMP, FLECH(UTC), ISOPUNSH (SOLLAC).
Cette méthode fournit des résultats approchéspar rapport aux modèles complets, mais elle
Les simulations
resteun outil très rapide et très utile dans la phasedu pré-dimensionnement.
'complets' correspondentà un schémade résolutionincrémentale.
ou les modèlesqualifiés de
L'approche incrémentale,par opposition à I'approche simplifiée, permet un suivi précis de
l'évolution des différents phénomènesrencontréspendantles différentesphasesdu processus.
Les approches incrémentales peuvent être basée sur une méthode statique implicite
(ABAQUS, MARC...), statique explicite ou dynamique explicite (OPTRIS, PAM-
STAMP,RADIOS, ABAQUS...).En raison de la complexitédes phénomènesà simuler
(couplagede plasticité,contact...),le temps de calcul resteun facteur très pénalisantdans
cesapproches.
3. Modélisations mécaniquesdu procédé de I'emboutissage
La simulation numérique de I'emboutissagenécessiteune description très fine des

phénomènesphysiqueset mécaniquesrencontréspendantles différentesphasesdu procédé.
L'introduction des lois de comportementdu matériau,aussiréalistesque possibleset prenant
en compte des processusphysiquessusceptiblesd'intervenir, est indispensable[Lemaitre/85]
[Gelir/90]. Dans la littérature de la simulation des procédésde mise en forme, differents

modèles ont été adoptéspour décrire le comportementdu matériau à emboutir. Il est difficile
de classerces approchespar méthodologie car la mise en forme des métaux fait intervenir des
problèmesdivers nécessitantdes approchesvariées.Néanmoins,on va citer danscettesection
cellesqui sont les plus utiliséesdansles simulationsnumériquesdu procédé.
3.1. Modèle riside plastique
Dans ce modèle, le comportement est supposé être parfaitement plastique; les

déformation élastiquessont négligées.On peut citer le modèle de [Lee et KobayachilT3)
[Kobayachi et Kim/78] [Wagoner et all88 ]. Dans d'autre modèlesI'aspect d'anisotropiea

été développé[Gotoh et Ishise/78] [Gotoh/80]. Le problème de contactet de frottementest
-15-
Chaoitre I Etude bibliographique
pris en compte, au début, dans le modèle rigide plastique de Zienkiewicz et al lZienkiewicz et

all78l. Dans un modèle rigide plastique, les équationsconstitutivesimpliquées sont plus
simples, comparéesaux autres approches.La formulation élémentsfinis conespondanteest
une formulation généralementen vitesse et les contraintes sont liées directement aux taux de
déformationset ils ne sont pas réactualiséspar leurs incrémentsmais calculés à la fin de
chaquepas. L'approximation faite dansce modèle peut êtrejustifiée dans le cas de forgeage
où le niveau des déformations plastiques atteintes est très êlevé, elles peuvent atteindre des
centainespour centscomparéesaux déformationsélastiques.Dans I'emboutissageà froid, les
déformationsélastiquesne sont pas négligeables,ainsi que les phénomènes
du retour élastique
et des contraintesrésiduellesne peuvent être accessiblesfacilementpar un tel modèle. Par
contre,avecdestempsde calcul de2 à 3 fois moins élevés(à discrétisations
égales),elles sont
particulièrement efficaces pour la simulation du forgeage où le taux de corroyage et le
remplissagedes formes en matrice sont principalementrecherchés[Brunetlg}].
3.2. Modèle viscoplastique
En général, le comportement viscoplastique est modélisé par une loi de Norton-Hoff
[Germain/85] [Lochegnieset all89] [Brunet/90] [Chenot et all92l [Chenot/93] [Hacquid96]

[Chenot/97] lRachik et Roelandt/98][Brunelot/99].Ce modèle est généralementutilisé dans
les procédésde mise en forme à chaud ou pour des matériaux dont I'aspect de viscosité est
dominant comme le cas des polymèrespar exemple.Les remarquesà propos de ce modèle
sont semblablesà cellesfaitessur le modèlerigide plastique
3.3Modèlesuperplastique
Certains alliages possèdent ce type de lois de comportements.Ces alliages

perrnettentdes élongations qui varient entre 200 et 1000% et peuvent atteindre, dans
Parmi les simulationsbaséessur ce modèle,on peut citer lZhang
certainscas,5000%o.
et all861 [Wood et all89\ [Bellet et all9Ùl. La superplasticitédu métal, dans une
certaine gamme de températureset de vitesses de déformation, permet de déformer,
sous des contraintesfaibles (pressiongazeuse),le métal et lui faire épouserla forme
-16-
Chapitre I Etude bibliographique
d'un moule. Ce type de loi de comportementest adoptésouventdans les simulations

desprocédésde mise en forme par gonflage.
3.4 Modèle élastoplastique
Dans la littérature, plusieurs approches sont basées sur ce modèle. Les premières
simulationsbaséessur ce modèle étaientlimitées à la théorie des petitesdéformations,et les
conditions de contact n'étaient pas prises en compte [Wifi/76] fWang et BudianskylTSl
[Duroux et Batoz/9}]. Ensuite, des modèles prenant en compte l'aspect de grandes

ont été développés[Tang/81] [Brunet/90][Gelin et al/95].Dans
déformationsélastoplastiques
le cas de I'emboutissage,un modèle élastoplastiquesemblele plus réalistepour modéliserle
comportementdu matériau. D'autant plus qu'il permet un accès facile aux phénomènesdes
contraintesrésiduelleset du retour élastique.La formulation d'élémentsfinis correspondante
à ce modèle est une formulation en déplacement,le plus souvent, lagrangienneréactualisée.
En élastoplasticité,la difficulté principale vient de la nécessitéd'une configuration de
référenceet du grand nombre des variables d'état impliquées dans la formulation. Par ailleurs,
la loi de comportement,sous sa forme incrémentale,doit être intégréeet transforméedans le
principe variationnel.L'approche élastoplastiquesemble la plus utilisée pour la simulation
des procédés de mise en forme des produits plats tel est le cas pour I'emboutissage
IBrunet/90].
4. Elémentsfinis utilisésdans les simulationsnumériquesde I'emboutissage
Les modèles de simulation numériquesde I'emboutissagedépendentfortement du tlpe

des élémentsfinis utiliséspour discrétiserla tôle. Dans la littérature,on trouve des approches
membranes,approches2D, approchessolideset approchescoques.
4.1 Approche membrane
Dans cette approche,seuls les effets des composantesdes efforts dans le plan de la tôle
sont pris en compte.Plusieursmodèlesbaséssur cette approcheont été proposés[Kobayachi
et Kim/781 [Wang et Budiansky 178] [Wane/8a] [Wood ell89] [Duroux et Batozl90f
[Brunet/9O]fBellet et all90] [Guo er all92]. Ces approchesont permis aux auteursde prévoir
-17 -
Chapitre I Etude bibliographiaue
la répartition des déformations et des contraintes ainsi qu'une bonne estimation des efforts
exercéssur les outils avec des temps de calcul et des espacesmémoire faibles par rapport aux
autresapproches.Mais le problèmequi se poseest qu'une approchepareillene permetpas de
prévoir les phénomènesdes plis et les phénomènesde cisaillement qui sont très importants en
emboutissage.
4.2 Approches bidimensionnelleset approchessolides
Dans le but de simplifier la simulation numériquede I'emboutissage,plusieursapproches

se sont limitées à une étudebidimensionnelleou axisymétriquedu procédé,ce qui a permis
d'aboutir souvent à de bons résultats.Dans toutes ces approches,les formes des emboutis
étudiées étaient simples. Ce qui n'est pas toujours le cas, surtout lorsque on parle de
I'industrie automobile.Les approchestridimensionnellesbaséessur des élémentsvolumiques
(solides)sont plutôt utiliséesdans la simulation des procédésde mise en forme des produits
massifstels que le laminageet la couléecontinueou le forgeage.Les modèlesbaséssur cette
dernière approche ne sont pas nombreux. Ceci s'explique par les temps de calcul très
importants généréspar ces types d'éléments.Par ailleurs, vu la faible épaisseurdes tôles
utiliséesdans I'industrie de I'emboutissage,une approchecoque semblela plus réalistepour
une simulationde ce procédé.
4.3Approchecoque
Dans cette approche,les effets de flexion dansle plan de la tôle sont superposés
aux effet membranes. Plusieurs simulations basées sur cette approche ont été
développées[Boivin/9O] [Onâte et Saracibarl92][Gelin et all95l [Picart et all99]. La
formulation coque semble la formulation la plus efficace dans la discrétisationpar
élémentsfinis destôles destinéespour la mise en forme.
-18-
Chapitre I Etude bibliosraPhique
5. Problème de contact et de frottement
fnfioduction
Les problèmesdes procédésde mise en forme incluent souventdes conditionslimites

de tlpe contact/frottement [Germain/85] [Durvillel93] lHacquin/961 lChenotl9Tl fFourment
et all97l. En effet dans les opérationsde mise en forme, les matériaux se déplacentpar rapport
aux outils le long du processus,ce qui génèredes efforts de pression et de frottement entre le
produit et les outils. Pour cette raison, la prise en compte des conditions de contact est
indispensabledans la simulationnumériquede ces procédés.Du point de vue numérique,la
simulations des procédésconstitueune tache délicate en raison du couplagede différentes
non-linéarités(géométrie,plasticité,contact).
Le problème de contact reste le problème le plus sensible à cause des nombreux
conditions à satisfaire et I'attention à accorderpour assurerla convergencede la solution avec
un temps de calcul raisonnable surtout en 3D où l'augmentation du nombre des points de
contact génère une réduction considérablede la taille des incréments de temps fNagtaal et
Taylor/9l1
Dans la simulation numérique de ces procédés, le problème de contact est
généralement assimilé à un contact entre solides déformables qui représente la matière à
déformer, et des solides indéformablesou obstaclesqui représententles outils (poinçon,
matrice, serre flan.)
Pour la résolution du problème de contactdivers types de méthodessont utilisés. On

peut citer la méthode de Lemke [Klarbring/86] [Pinto/9O], la méthode de pénalisation
[Jear/93] [Cumier et AlarV8S] [Brunet/88] [Ben Dhia/89] [Lochegnieset all89l,la méthode

des multiplicateurs de Lagrange[Bathe et Chaudhary/85][Bath et Chaudhary/86][Courtois et
Ricard/881[Combescureet all88] et la méthodedu Lagrangienaugmenté[Ju et Taylor/88]
fFrachome et all89l [Simo et Laursen/92]. Une autre classede méthode consiste à introduire
des éléments finis spécifiquesintégrants toutes les conditions de contact dans leur définition
fGrolad et Critescu/S3] [Charlier et Cescotto/88].Toutes ces méthodesprésententleurs

avantageset aussi leurs inconvénients.Les schémasde ces différentes méthodes remplissent
à la prise en comptedesphénomènesde contactet
tous à leurs manièresles étapesnécessaires
de frottement[Jeanet Touzot/88],à savoir :
-19-
Chapitre t Etude bibliosraphique
- La détectiondu contact
- La détectionde la perte de contact
- Le respectde la condition de non interpénétration
- La détectiondu glissementou de I'adhérence.
5.1. Positiondu problème
On considèreun corps S occupantle domaine f) et ayantpour frontière ôC)

ôf)=ôÇùr+ôÇ20+ôQ"
Le solide S est soumis à une densitéde force surfacique f sur la frontière ôÇ)o et une densité
de force de volume / dans Q . La frontière ôf), est supposéeêtre à déplacementimposé z,
(figure 1.2).
Figure (1.2) Corps déformable en contact avec une surface rigide
Sur la frontière ôO., le corps est susceptibled'être en contactunilatéral.Le corps S serait
alors soumis à une densité surfaciquede forces due à la réaction de la surface igide f ,.
Généralement,l'étude mécanique d'un tel problème revient à rechercherle champ de

déplacementcinématiquementadmissible U du solide engendrépar le chargementdécrit ci-
dessus.
Dansun cadrestatique,le problèmerevient à résoudrele systèmed'équationssuivant :
-20-
Chaoitre I Etude bibliograohiaue
divo+T=O sur ç) (1.1)
o.n=l sur ôC)F (r.2)

o.n=f; sur ôC). (1.3)
u=ud sur âÇ), ( 1 .4)
A ces quatre conditions s'ajoutent la loi de comportementdu matériau qui peut être élastique,
plastique,élastoplastique
ou autre.
5.2. Conditionsde contact
Le problème de contact génère d'autres conditions à prendre en compte par rapport

aux autres problèmes mécaniques.Les conditions de contact, dans le cas d'un contact avec
frottement,portent simultanémentsur le déplacementet les forces de contactinduites.
Corps déformable
Corpsrigide
Figure (l.3) conditions de contact unilatéral
Les points (ou les næuds après discrétisation par éléments finis) appartenantà la surface de
contactdoivent satisfaireaux conditions suivantes:
* Conditiondenon-interpénétration
desdeuxcorps:
h>0 (1.5)
* Conditionde contactunilatéralou conditionde non adhésion:
- 21 -
4.È0 (1.6)
L'effort a pour effet de comprimer le corps au point de contact.

* Condition de complémentarité
F,lt=o (r.7)
Lorsqu'il y a contact,l'écarth estnul et une réactionde contactest générée.Pour un écartnon
nul, la réactionestréduiteà zéro.
5.3. Lois de frottement utiliséesdans les simulations numériques
Dans les procédésde mise en forme à froid, ce qui est le cas de l'emboutissage,le
phénomènede frottement joue un rôle important. Plusieurs simulations du procédé ont tenu
comptede ce phénomèneen proposantdesprocéduresde gestionplus ou moins complexes.
Dans la pratique, il existe trois catégoriesprincipales du frottement :
o frottement quasi statique

o frottement hydrodynamique
r frottement avecpénétrationet usure
Dans I'industrie desprocédésde mise en forme, le phénomènele plus connu fait partie plutôt
de la premièrecatégorie.
5.3.a. Modèle de Coulomb
Ce modèle met en jeu un coefficient de frottement qui est le rapport de l'effort

tangentielà I'effort normal exercépar un solidesur l'autre :
F,<ÀF,I
si F,</,r,1 ona v,-0 (1.8)
si F,=/,F,1 rI = 1> 0 tel que v, - - lF ,
11
Chapitre I Etude biblioqraPhique
p est le coefficient du frottement, v, est la vitessetangentielle du point considéré, F, et F,
représententrespectivementla composantetangentielle et la composantenormale de la force

du contactau point étudié.
Figure (1.4) loi dufrottement : Modèle de Coulomb
Cetteloi peut être écriteaussien fonction de la vitessetangentiellerelative(ni"')
F,=-
lrll,lffi si v',"t+O (t.e)
F,.tlF,l si vi"/=g (1.10)
On remarqueque cette loi présenteune discontinuité, ce qui va créer des difficultés lors de la
simulation numérique. Pour surmonter ces difficultés, des techniquesde régularisation sont
souventemployées:
figure (1.5); régularisation de la loi defrottement.
-23-
Dans les industriesd'emboutissage,on utilise rarementles lubrifiants: ceci peut justifier une
hypothèsed'invariance du coefficient de frottement.
5.3.b. Modèle de Tresca
Cetteloi s'écrit sousla forme linéairesuivante:
si lF,l.#Po
' , vl,=o (1.11)
a/J
Si lF,l=+po y ,- - t F ,
it= 1>0 telque (r.r2)
{J
m représentele coefficient de frottement de Tresca.
lnl Coulomb
lF,l
Figure (1.6) Modèles defrottement: Coulomb, Tresca.
On remarqueque le modèle de Trescaprésenteun seuil de frottementconstant.Les forcesde

contacttangentiellesinduitespar le frottementsontindépendantes
de la réactionnormale.
5.3.c. Modèle de frottement non local
Ce modèle a été proposépar J.T Oden et E.B Pires fOden et Pires/83][Oden et Pires
/88] [Oder/88]. Au contrairedes lois de Trescaet de Coulomb , qui déterminentla valeur de
la force tangentielle{ à partir de I'effort normal de contact { en un point P, dansce modèle,
-24-
elle est déterminéeau voisinagede ce point. Les auteursont proposéaussiune loi non locale
et non linéaire par rapport au déplacementtangentiel u,.
Vu la complexité de la formulation, on se limite à cette brève présentation.Par

ailleurs, ce modèle est rarement utilisé dans la gestion de problème de contact dans les
simulationsnumériquesdesprocédés.
La mise en ceuvre du modèle de Tresca reste aisé du fait de la simplicité de son

écriture en comparaison à la loi de frottement de Coulomb, qui est plus élaborée dans sa
définition.
La loi de Coulomb décrit bien le phénomènede frottement si ln,l n'est pas très
important, par contre la loi de Tresca est surtout proposé dans le cas otr lF, I est grand. Pour
les procédésd'emboutissageles pressiottr lF,l de contact sont beaucoupmoins importantes
que lors du forgeage et du filage. Par ailleurs, en vertu des essais industriels de
I'emboutissage,la loi de Coulomb est mieux adaptéepour représenterle frottement entre
métauxque la loi de Tresca[Avitzuri83].
5.4. Méthodes de traitement des problèmes de contact et de frottement
La résolutiond'un problèmede calcul de structuresnon-linéairesanstenir comptedes

conditionsde contactrevient à la recherchede la solutiond'un systèmede tlpe:
lrorl tuFlrl ( 1 .l 3 )
où K est la matricede rigidité qui dépenddu vecteurde déplacementsolution du problèmeU.

Le vecteur F, est un vecteur secondmembre qui représenteles forces extérieures.En tenant
comptedes conditionsde contact,la dernièreéquationdevient :
Ktul tuFtFl.tRl (1.14)
Le vecteurR représenteles réactionsde contact.

Ces réactions de contact sont indéterminées a priori et constituent des inconnues
supplémentairesau problème à résoudre.
Plusieursméthodesont été proposéespour traiter le problèmede contact:
-25-
5.4.a Méthode de Lemke:
La méthode de Lemke [Klarbring/86] fRaous et all88l [Pinto/90] [Chertier/97) traite

en premier la zone de contact ( condensationsur la surface de contact) par opposition aux
autresméthodes qui traitent globalement I'ensemble des corps en contact et leurs interfaces.
La résolution du problème de contact,à I'aide de cette technique,revient à :
Trouver w€RP , zeRP tels que :
w_Mz=q (1.15)
w; >0 (1.16)
zi>O i=|,....,p (1.17)
wizi=0 (l'18)
où q est un vecteur donné de Rp et M une matrice carrée donnée et de rang p. Il est
nécessairede passerpar une techniquede condensation[Lebon/88]du système:
rxrurl
[l:]=r{t] (r.1e)
où les indices / et c indiquentrespectivementles næudslibres et les næudsen contact.Ce qui

revient à résoudre:
k.ol [u"]-h,F[r.] ( 1. 20)
de K et F.
K" et F. sont les formescondensées
Par identificationavec(1.15),(1.16),(1.17) et (1.18),on a :
W=-R")Q , z=-\")0, q=F- , M=K'
La méthode de Lemke est une méthode mixte puisqu'elle comporte la résolution en

déplacementet en force de contact. Cette méthode présenteune simplicité et une performance
en vitesseet en qualité de résultats.Néanmoins,pour des problèmesoù le nombre de points
de contact est grand (>200) il peut apparaîtredes phénomènesd'accumulations d'erreurs dus
principalementà la condensation[Kallala/93].
-26-
Chapitre I Etude bibliosraphiaue
5.4.b.Méthode des multiplicateurs de Lasranse
Cette méthode est utilisée dans les problèmesd'optimisation sous contraintes.Elle

consiste à introduire des paramètressupplémentairesdits multiplicateurs de Lagrange. Elle
revient, en fait, à introduire des contraintessupplémentairessur les næuds en contact dans le
systèmematériel global. Plusieurs algorithmes basés sur cette méthode ont été développés
[Bathe et ChaudharylS1][Bathe et Chaudhary/86] fBathe et Mijailovich/88] fCourtois et

Ricard/881[Combescureet all88].
La condition cinématiquede contactpeut s'écrire en un point de contact sousla forme
suivante :
p=Q".u"+6=0 (Lzr)
Où â est la pénétrationinitiale du point de contact, p est la pénétrationactuelleliée
au déplacementu"du point de contactpar I'intermédiaired'un opérateurQ, lzhon!931. Pour
un ensemblede N points de contactl'équation (1.21) s'écrit d'une manièreplus généralesous

la forme :
p=eU+A (r.22)
Dans la méthode de multiplicateurs de Lagrange,la fonctionnelle à extrémaliserest remplacée
par:
n, e,^)+u'KU -u'F+^1(QU
+^) (r.23)
où ^ est un vecteur de même dimension que P. Ses composantess'appellent les

multiplicateurs de Lagrange. Le problème d'extremum revient donc à trouver U et Â tels
que:
(r.24)
#=o
ft=o (r.2s)
L' équation(I.24) donne:

KU-F+QÂ=0 (r.26)
et l'équation(1.23)donne:
QU+Â=0 (r.27)
En combinantles équations(1.26) et (I.27), on obtient :
-27
Chapitre I Etude biblioqraohique
KrUr=4, (r.28)
K,=[6]'] ,q=(-ï),u,=(X)
avec
En résolvantl'équation (1.26), on obtient le déplacementU et le multiplicateurde Lagrange
Â. On peut constaterque dans cette méthode,les multiplicateursde Lagrangeont le sens
physique des forces de contact. La méthode des multiplicateurs de Lagrangeest une méthode
mixte puisqu'ellecomportela résolutionen déplacementet en mutiplicateurs.
Comme on peut le constater,I'utilisation desmultiplicateursde Lagrange,multiplie le

nombre d'inconnuesdu problème. L'apparition d'un bloc de zéros dans I'opérateurtangent
pose des problèmes numériques dans la résolution et exige I'utilisation des algorithmes
particulierstel que la méthodede type gradient[Chenaiset Rousselet/88].L'avantagede cette
méthode est que les conditions aux limites sont exactementsatisfaiteslZhon!93|.
5.4.c. Méthode de pénalisstion
La méthode de pénalité est parmi les méthodesqui sont courammentutilisées dans les codes
d'élémentsfinis pour sa grandefacilité de programmation[Hunek/93]. L'étude du problème
de contactpar cetteméthodepermet de passerd'un problèmed'optimisation souscontraintes
à un problèmed'optimisationsanscontrainte.
Cette méthode consiste à introduire, dans la fonctionnelle du problème, la fonction de
pénalisationsuivante:
r o=lP'aP (t.zg)
Où a est une matrice diagonale de composantes au qui constituent les coefficients de
pénalisation.La fonctionnelle à minimiser devient donc:
fl o=ll+ x o
:lu'KU -u 'r +!Y'aY (1.30)

Et par suite, le champ de déplacementU solution du problème est le champ qui minimise la
fonctionnelle fI o :
-28-
ôfl _
-ffi=u (r.31)
En substituant(1.22) et (1.30)dans(1.31),on obtient l'équationsuivante:

KoU=Fo 0.32)
AvecKp=K+Q'aQ et F,=F+Q'aA
En résolvant l'équation (1.32) on obtient le déplacementU. Les forces de contact sont

donnéespar la relation :
F,=d P (U ) (1.33)
Cette équationest résolue,après,au moyen d'une méthodede Newton-Raphson.Ce

qui permet de déterminer l'état des næuds de la zone de contact à chaque itération. Dans la
méthodede pénalisationle déplacementest la seuleinconnue.
En ce qui concerneles coefficientsde pénalisation,ils doivent être assezgrandspour
minimiser l'interpénétration, tout en gardantun bon conditionnementde la matrice de rigidité.
Ceci n'est pas toujours évident, en fait les conditionsde contactne sont satisfaitesque pour
des valeurs infinies des coefficients de pénalisations ce qui engendre un mauvais
conditionnement de la matrice de rigidité. Cette méthode entraîne une réduction du pas de
temps pour un schémad'intégration explicite. Par ailleurs, elle ne permet pas de vérifier
exactementles conditionsde contactet dépendfortementdu choix du coefficientde pénalité :
un mauvais choix entraînedesrésultatsphysiquementinacceptables.
5.4.d. Méthode du Lasraneien ausmenté
Cette méthode est une combinaison entre la méthode de multiplicateur de Lagrange et la

méthodede pénalisation.Elle consisteà satisfaireles conditionsd'équilibre tout en gardantun
bon conditionnementde la matrice de rigidité et ceci en choisissantdes coefficients de
pénalisationfinis ou moins 'fortement' que pour la pénalisationseule.Elle présenteun autre
avantage: il n' y a pas une augmentation des paramètresétant donné qu'elle revient à une
pénalisationdu multiplicateurde Lagrangequi n'est plus considérécomme inconnu lors de la
résolutionmais qui seraaugmentéau coursd'un processusitératif.
Dans cetteméthodel'expressionde l'énergiepotentielletotale devient:
1u=fI ,+!t{etr ( 1.34)
-29-
Chapitre I Etude bibliogranhiaue
Où a estune matrice diagonalequi contientles paramètresde pénalisation€,,, i=l,N
La stationnaritéde fI ,, permet d'écrire :
ôfl u =o (1.35)
âU
ôfIu =o ( 1.36)
aL
Les équation(1.34) et (1.35) donnentl'équationdu problèmeà résoudre:
KrrU r='F, (1.37)
avec
K,"=[ô:'1
Le multiplicateur de Lagrange est réactualisé après chaque pas d'équilibre (algorithme

d'Uzawa) ou aprèschaqueitérationd'équilibre [Simo et Laurser/92].
Dans le deuxième cas, il n'y a pas de boucle, mais la réactualisationdu multiplicateur de
Lagrange modifie la convergence quadratique du schéma consistant de Newton-Raphson
fSimo Laursenl92f.
Dans le contexte du problème de contact avec frottement une formulation à trois champs
indépendants( Hu Washizu ) a été proposéepar Papadopouloset Taylor [Papadopouloset
Taylorlg3l. Dans cette formulation, un contact entre éléments est postulé au lieu du contact
entre les næuds. Ce qui introduit une fonction considéréecomme variable indépendantedans
la formulation.
5.4.e.Méthode des élémentsd'interface
Cette méthodeest baséesur I'introduction des élémentsaux sensélémentsfinis entre

les zones susceptiblesd'entrer en contact fGrolad et Critescui83] [Charlier et Cescotto/88].
Les éléments de contact appartiennent à la fois aux deux solides en contact et ont un
comportementprédéfini. Ces élémentsd'interfaceréagissentau signede la contraintenormale
moyenne. Dans le cas d'une compression(valeur négative de la contrainte), le rôle de
l'élément est maintenu,sinon il y'a rupturede contact(traction).
Dans le casde contactglissant,la matriced'élasticitéde l'élémentde contacts'écrit :
-30-
Iooo
hFl oEo
1000
En présencede contactavecfrottementla matriced'élasticitédevient:
..1-ooo'l
hFl et défini parr,n=Gy,n
0 E O Ioù G estun termede frottement
L0oc_l
Au niveau de la matrice de rigidité, la prise en compte des éléments de contact se
traduit par I'addition d'un bloc supplémentaire.
Ainsi, la matrice de rigidité K de la structure
aura la forme :
t"t=FLIJ
Kn et Ke étant les matrices de rigidité de deux corps A et B susceptiblesd'entrer en contact et
Kc la matricede rigidité desélémentsde contact.
Cette approche de contact donne des résultats satisfaisantspour des géométries de
zonesde contactsimples.Elle présentedes contraintesau niveau de la discrétisationen espace
et son application est très limitée. Cette méthode a été implémentéedans plusieurs codes
industrielstels que CASTOR [Grolad et Critescu/83]et ABAQUS [ABAQUS/5.3].
6. Méthodes des résolutions numériques
Introduction
L'établissementd'une méthode de simulation numérique des procédésde mise en

forme se fait sur trois étapes:une formulation théorique du problème qui permet d'aboutir aux
conditionsmécaniques(équilibre, loi de comportement,contact,...)puis une discrétisationdu
problème (exemple méthode des éléments finis), ensuite une résolution des systèmes
d'équations non-linéairesobtenus. Ces demières années,la méthode des élémentsfinis a
connu un essor considérable,et est maintenant utilisée dans la plupart des codes et les
méthodes numériques pour les sciences de I'ingénieur, elle permet d'effectuer une
discrétisationspatialedu problème à résoudre.De nombreux ouvragesexistent sur le sujet:
-31-
fDhatt et Touzot/84] [Batoz et Dhatt/90] [Imbert/9l] [Zienkiewicz et Taylor/9l] [Bathe/96].

Malgré le nombre important des codes de simulations numériques existants, le nombre des
méthodes et des techniquesde résolutions numériques des équationsnon-linéaires reste très
limité.
6.1 Méthodesclassiques
La techniquela plus connue est celle de Newton Raphsonou techniquede prédiction-

correction. Cette technique procède en deux étapes; la première consiste à linéariser le
problème non-linéaire de départ et donc à prédire une nouvelle solution, tandis que la
deuxième consiste à corriger le résidu d'équilibre global de la structure par itérations
successives.Ces méthodessont bien adaptéespour la résolutionde toutes les non-linéarités
possibles: plasticité, contact,.... En revanche, elles n'exploitent pas les particularitésde
chaque type de non-linéarité, et le choix des paramètresde contrôle demandeun utilisateur
expérimenté.Le choix de la taille des pas est très délicat dansces techniques; un incrément
de temps très grand n'assurepas la convergencede la solution, et un incrément très petit
multiplie remarquablementle coût du temps de calcul. Cependant,plusieurs articles ont été
consacrésà I'amélioration de la convergencede ces algorithmeset à I'adaptationdu pas du
calcul fCriesfield/81] [Criesfield/83].
Depuis quelquesannées,quelqueséquipesde recherchetravaillent,non pas seulement
sur I'amélioration de ces algorithmes, mais aussi sur le développementde nouvelles
techniqueset algorithmes de résolution des problèmes non-linéaires,en vue de proposer des
outils précis, robusteset fiables. L'équipe de Ladevèzetravaille sur une méthodede calcul
non-incrémentaledite 'méthode à grands incréments de temps lLadevèze/85] [Ladevèzel9l]
lLadevèzel96l
Malgré I'augmentationrapide des puissancesdes calculs informatiques,les logiciels
de simulation de I'emboutissagenécessitentencore des temps de mise en données et de
calcul très importants et ne permettentdonc, au bureaudes méthodes,que de procéderaux
ultimes vérificationset validationsavantla fabricationdesoutils [Lachatet Ravier/99].
6.2. Les méthodes asvmptotiques numériques
La Méthode asymptotique numérique (M.A.N.) est une famille d'algorithmes qui permet de
résoudredes problèmes non linéaires grâce à I'associationdes techniquesde perturbation et des
-32-
ChapitreI Etudebibliographique
méthodesnumériques telles que la méthode des élémentsfinis. En effet, la M.A.N. permet de chercher
les branches solutions sous la forme de séries entières par rapport à un paramètre "a" que nous
préciseronsplus tard. De cette manière,nous transformons le problème non linéaire de départ en une
suite récurrentede problèmes linéaires,bien posés,admettanttous le même opérateurtangent. Ensuite,
ces problèmes linéaires sont résolus numériquement,généralementpar la méthode des éléments finis.
Comme le rayon de convergence des séries est fini, une procédure de continuation a été proposée
[Cochelin/94] permettant de définir la longueur de chaque pas asymptotique de sorte que le résidu
reste très faible. Chaque fin de pas servira alors comme point de départ pour le pas suivant. La
longueurde pas est par conséquentdéterminéea posteriori; ceci n'estpastoujoursle cas desméthodes
itératives habituelles. Ainsi, nous avons une longueur de pas adaptativeet déterminéede façon simple
et optimale : c'est le point clé de I'efficacité et de la robustessede la méthodeasymptotiquenumérique.
L'idée est donc de tracer numériquementla branche solution non pas par une séquencede points mais
par une séquencede morceaux.
En partant de la formulation en déplacement d'un problème tridimensionnel d'élasticité
en non linéaire géométrique avec une loi de comportement linéaire, l'énergie potentielle est
stationnaire pour tout champ de déplacement cinématiquement admissible et s'écrit sous la
forme suivante:
ôII = I'y|ru):D:fu
(u)dv- )"Pe( ôu)=g (1.38)
Où D est la matrice de comportement, u le champ des déplacements,)" w paramètre de

chargementet Pe est le travail des efforts extérieurs.y est la déformation de Green-Lagrange
que nous pouvons décomposeren une partie linéaire et une autrepartie non linéaire:
(u)+y" çu1=NV
74u1=yt u+'Vu)+N'v uv u;
(1.3e)
6r G)=y'(6u )+ 2 y'' (u,6t)
Le principe de base de la technique de perturbation consiste à chercher les branches de
solution du problème (1.38) en développant u et 2 en séries entièrespar rapport à un
paramètrede cheminbien choisi " a". Ce paramètre,représenteune inconnuesupplémentaire,
il peut être défini de la même manière que les paramètresde contrôle pour les algorithmes
itératifs classiques(chargeimposée,déplacementimposé ou longueurd'arc imposée).Ainsi,
nous pouvons écrire :
-JJ-
P p
Sat
u(al=>4 oui, ).@)=l a').i (1.40)
i=0 i=0
En injectant(1.39) dans l'équation (1.38) et en identifiant les termes obtenus suivant la

puissancede " a ", nous obtenons une suite récurrente de problèmes linéaires que nous
pouvons écrire à I'ordre p quelconquesousla forme suivante :
Lt(uo)=)rpaPit (1.4r)
où Lt représenteun opérateurtangentde rigidité qui dépenduniquementde la solution initiale
(lo ,4), F le vecteur des forces appliquées,Ei un vecteur secondmembre ou un vecteur
résidu qui contient des termesdépendantdes ordres inférieurs à p .
(1,(u ,),6u)= I't y{u p)+2y,r(uo,up):D:y(6u)dv +
* I ' ( r , , ( uo ) +2 y , , ( u e , u
o ) : D: T , t ( uo . 6 u ) d v (r.42)
* I ' ( y ( u p ) + 2y a ( u o , up ) : D: y n ( u o , 6 t ) d v
(tr',âu)=P"(âu) (r.42)
-
)= I E ,I r , , ' , ) * f o ' r , 1 ( u . , u , - , ) :)2
(L''0,6u ,vn, ( n , - , , 6 u ) d v
_JL'y ,,(, ,,, o-):D:(T,(6u)+2y,,(uo,6t))dv

p-l
r=l (1.43)
Remarquonsque I'expression du vecteur plt îait intervenir une double sommation due à la
formulation de départ qui est cubique. Cela peut avoir une influence importante sur le temps
de calcul.
Après avoir établi la condition d'équilibre à I'ordre p, nous procédons à la
discrétisationdu problème pour le résoudrepar une méthode des éléments finis de type
déplacement.Le problèmeglobal discrétisés'écritsousla forme suivante:
-34-
t K r l { uo YA , { "} . b ; ' } (r.44)
où [K4 désignela matrice de rigidité tangenteécrite au point représentantla solution initiale,
{uo} tr vecteurdéplacementdiscrétiséà I'ordre p .

Le but de cette méthode est donc de parcourir une grandepartie de la branche solution
en n'inversantqu'uneseulematrice [Kr] communeà tous les problèmeslinéaires(1.8), et de
reporter tous les termes non linéaires sur le vecteur tri'l qui dépenddes ordres précédant
I'ordrep.
Parmi les premiers travaux les plus cités dans la littérature sur I'application des
méthodes de perturbation,on trouve celui de Signorini [Signorinii30] qui a explicité les
relations entre l'élasticiténon linéaire et l'élasticitélinéariséeen écrivant le problème non
linéaire sous forme d'une suite de problèmesd'élasticitélinéaire. Vient ensuiteKoiter qui a
écrit en 1945 sa thèse sur la théorie des bifurcationsen élasticiténon linéaire dans laquelle,
pour caractériser la sensibilité aux imperfections des structures élastiques,il a utilisé une
méthode de perturbation tronquée à I'ordre 2lKoiterl45]. Ces travaux ont été repris plus tard
par Budiansky (1974) [Budiansky/74]puis par Potier-Ferry[Potier-Ferry/87].A noter que la
majorité des travaux sur les méthodes de perturbation pour le calcul de structure trouve ses
origines dans la théorie de Koiter. Jusqu'ici, les méthodesde perturbation, de part leur aspect
analytique,ne résolvaientque desproblèmesà géométriesimple.
Thompson et Walker fThompson et walker/68], ont été les premiers à associerces

techniquesà la méthode des éléments finis pour l'étude des problèmesde mécaniquedes
structures.Ils ont élaboréune théorie généralepour résoudreun problème non linéaire discret
par une méthode de perturbation qu'ils ont appliqué ensuiteau problème d'équilibre de poutre
élastique. Ils ont obtenu alors une successionde problèmes linéaires par dérivations
successivesdes équationsde départ par rapport au paramètrede développement<<c>>.Ils ont
ensuite déterminé les sept premiers termes de la série pour une poutre en flexion avec une
non-linéaritégéométrique,et ils ont comparéles résultatsobtenusavec la solution analytique
exacte. Ils ont montré que plus I'ordre de troncature augmente, meilleure est la précision
asymptotique. En revanche, les auteurs ont conclu que ces méthodesne peuvent traiter que
-35-
'Walker
des problèmes avec des non-linéaritésmodérées.Plus tard, (1969) [Walker/69] a
publié un article dans lequel il a étudié la stabilité et le comportementd'un arc circulaire en
grands déplacements.Il a utilisé les méthodesde perturbation comme moyen de prédiction de
la solution qu'il corrigeait systématiquementà la fin de chaque pas par les procédures
itératives de Newton-Raphson.Pour ses développements,Walker est parti d'une formulation
cubique en déplacement,du même tlpe que (1.38) avec une troncaturede la série à un petit
ordre. Un certain nombre de travaux sur les applications de ces algorithmessont résuméspar
Gallagheren 1975 [Gallagher/75].Toutefois, I'auteura conclu que ces méthodesne sont pas
amenéesà se développerdansle cadredesélémentsfinis.
Les années80 ont connu une relance de I'intérêt des méthodesde perturbation. Ainsi,
Noor et ses collaborateursont associéces demièresà la méthodede réduction de Rayleigh-
Ritz fNoor et Peters/80][Noor/8l] fNoor et Peters/8l]. Partant du principe que la réponse
d'un problème en statique est décrite par un systèmed'équationsalgébriques,cette méthode
consiste à remplacer ces équationspar un systèmeréduit d'équationsavec beaucoup moins
d'inconnues.Ces auteursont utilisé la méthode de perturbation pour générerla base de Ritz.
Quelques avantages peuvent se dégager de ces nouveaux algorithmes, ptr exemple la

réduction de la taille du problème à résoudreet I'amélioration du rayon de convergencede la
solution asymptotique, perrnettant ainsi d'optimiser le nombre de décompositions de la
matrice de rigidité. Plusieursarticles montrent I'efficacité de cette méthode.Cependant,le
dans ces articles,alorsque le but recherchéest,
temps de calcul n'estpas discutésérieusement
justement,d'allier I'efficacitéà l'économiedu temps de calcul; voir I'articlede synthèsesur ces
méthodes,présentépar Riks en 1984 [Riks/8a].
En 1990,Damil et Potier-Ferry[Damil et Potier-Ferry/90]ont écrit le premier article

original qui traite des bifurcations perturbées à I'aide des méthodes asymptotiques
numériques. En effet, pour obtenir une efficacité numérique et une écriture simple du second
| ;l
membreF "' j jusq"'à des ordres élevés,ces auteursont proposéune formulation de départ
quadratique.Au lieu de partir de la formulation en déplacement(1.38) qui est cubique, Azrar

et al.lAzrar et all93f ont proposéde réduirele degréde non-linéaritéde u en introduisantle
champdes contraintescommevariableindépendante.Si nous considéronsl'élasticité3D, cette
formulation de base correspond à la fonctionnelle d'Hellinger-Reissner,qui est stationnaire
-36-
pour tout champ de déplacement cinématiquement admissible âu et un champ âS

quelconque:
el -'
u)*'ôS:(u)-' S:D :âSfv -ZP"(ôu)
(u,S)=J{s'a71
ôÉ1R (1.4s)
Les variables développéesen série sont donc le paramètrede charge 2 ainsi que le champ de
déplacement u et le champ des contraintes S que nous pouvons rassemblerdans un vecteur
mixte U=[X). L'avantagede cette écriture est que le problème à résoudreest écrit dans un
\ù/
cadre quadratique que nous pouvons mettre sous la forme suivante en fonction du vecteur
mixte:
R(U,2)=L (U )+Q (U'U )- \f.:0 (r.46)
où R est le vecteurrésidu, L un opérateurlinéaire et Q un opérateurquadratique

Nous avons les relations suivantes:
u)- D-':s)pv
(r tu ),âu)=Jf s'r,tau)+'âS:(,y,( (r.47)
(Ofu,u),ôu)=Itt :21,,(u,ôu :7,,(u,r

)+'âS )h ( 1.48)
Si nous appliquons maintenant la méthode de perturbation, nous obtenons à I'ordre p u n
problème linéaire qui s'écrit sousla forme suivante:
L , ( U ) = ) "o F + F i l (r.4e)
où L, correspondau même opérateurde rigidité explicité dans le cas de la formulation
cubique.Le nouveauvecteur non linéaire Ffl s'écrit quant à lui comme suit :
p
F i ' =-I qru,,up-,) ( 1. 50)
avec :
(Q{u,,uo-,),flJ)= o-,)Vu (l .51)

Jfs,'r,,{uo-,,ôu)+'So-.:y,,(u,,âr)+'ds:yn,(u,,u
a1
-Jt-
Contrairementà la formulation cubique (1.38) qui fait intervenirune double sommationdans
l'expressiondu secondmembre Fi' ,la formulation quadratiquen'a besoin que d'une simple
sommation.Le temps de calcul de {rJ'} ,.u. petit par rapport à celui nécessairepour le
calcul de la matrice de rigidité tangenteet par conséquentun pas asymptotiqueexige un temps

CPU du même ordre que celui correspondantà la méthodede Newton-Raphsonmodifiée.
Dans son travail de thèse, Zahrouni fZahrotnil98] a présentéune étude détaillée du
coût d'une formulation quadratiqueet son intérêt par rapport à une formulation cubique.
A I'ordre p, nous procédonsà une condensationdes contraintesen remplaçantla variablede
contraintes par la loi de comportement qui sera donc respectéeexactementen chaque point
d'intégration.Nous utiliserons par la suite une discrétisationhabituellepar élémentsfinis en
déplacement.
Cette idée fut vite généraliséeet son champ d'applicationélargi pour résoudredes problèmes
de coques en non linéaire géométriquedans le cadre des rotations modéréeslAzrar/931, mais
aussi pour le calcul des points de bifurcations [Boutyourl94] et la recherche des branches
post-critiquesde flambage[Boutyour et all93] [Vanucci et all98), ainsi que pour la résolution
deséquationsde Navier-StokesïnÆ6l [Ti et all96] [Cadou/97][Hadji/95]...
Comme la représentation en série est limitée par son rayon de convergence (propriété
intrinsèque), quelques algorithmes ont été développéspour améliorer la représentationde la
solution par série entière [Cochelin et all94] [Baker et all96]. Nous pouvons remplacer
I'approximation polynomiale par des fractions rationnelleset ainsi écrire :
u(a)-u(0)=2rt@)-o, (r.s2)
on ft@) sont des fractions rationnelles appeléesaussi approximantsde Padé.C'est une façon
simple d'améliorer le rayon de validité des séries;toutefois une fraction rationnelle présente
malheureusementdes pôles (des zéros au dénominateur) et des défauts (facteurs communs
entre le numérateuret le dénominateur).
Najah et al. [Najah et all98] ont présenté une revue critique des méthodes
asymptotiques numériques (techniques de bases réduites, représentation polynomiale et
approximants de Padé). Il en ressort q'une efficacité de ces algorithmes est obtenue avec un
-38-
ordre élevé de troncature des séries. Quant aux conclusions, elles prônent la représentation
polynomialeen raisonde sa fiabilité et de sa rapidité,ou les approximantsde Padéà condition
d'améliorer la fiabilité de cette technique.
Récemment,une procédure de continuation baséesur les approximantsde Padé a été
élaborée par Elhage-Husseinet al fElhage-Husseinet aV00] et implémentée dans de
nombreux algorithmesasymptotiques.
Les méthodes asymptotiques ont montré leur efficacité pour la résolution des
problèmes de mécanique des fluides [Cadou/97], plasticité [Braikat/95], viscoplasticité
[Brunelot/99], contact unilatéral [Elhage-Hussein/98]et coques en grandes rotations

[Zafuouni/97]. Récemment,les applicationsdes méthodesasymptotiquesont été élargi pour
traiter les problèmesde dynamiquefCochelinet Compain/991,problèmesde grandestailles,
et les procédésde mise en forme [Brunelot/99][Abichou et all}Il.
sous-structuration
-39
1.7.Conclusions
L'emboutissagecommeprocédéde mise en forme implique de nombreuxphénomènes

physiqueset mécaniques(plasticité,contact, frottement,...).La simulation numériquede ce
procédés'avèreindispensablepour prévoir I'emboutissabilitédesproduits et diminuer le coût
et les délaisde fabrication.Une telle simulation,pour qu'elle soit la plus réalistepossible,doit
décrire d'une façon très fine les différents phénomènesphysiques rencontrés le long du
processus.Lors de la simulation,ces phénomènesphysiquesse manifestentpar un couplage
de différentes non-linéarités: non-linéarité géométrique, matérielle, et non-linéarité de
contact.C'est ce couplagedes non-linéaritésqui rend la simulationdifficile à élaboreret qui
multiplie, par ailleurs, le temps de calcul des codes et des logiciels actuels.Les logiciels
destinés à la simulation numérique en mécanique utilisent souvent des méthodes
incrémentalesitératives pour la résolution des problèmes non linéaires. Ces méthodes
classiquesnécessitentdestemps de calcul très élevésquandil s'agit de problèmesde grandes
tailles et qui combinentde nombreusesnon-linéarités.Ce qui est le casdansla simulationdes
procédésde mise en forme en généralet de l'emboutissageen particulier.
Dans ce travail de thèse, nous allons présenter une simulation du procédé de

I'emboutissagebasésur les méthodesasymptotiquesnumériques,qui constituent,maintenant,
une alternativeintéressantepour la résolutiondes problèmesnon linéaires.Notre objectif est
d'alléger les temps de calcul qui présentent un grand problème dans les simulations
numériquesdu procédé.Nous testonsdonc, notre méthodedans la résolution des problèmes
qui combinent différentes non-linéarités.On note que dans ce travail, nous n'avons pas
I'intention d'élaborer un logiciel complet de simulation du procédé; ceci constitueraI'objet
de nos prochainstravaux.
-40-
Chapitre 2 Formulation du problème de contact entre un corps rigide et une coque
Chapitre 2
Formulation du problème de contact entre un corps rigide et une coque
fntoduction 42

2.1.1.Descriptioncinématiquede la coque 44
2.1.2.Equation d'équilibre 45
2.1.3.Loi de comportement 46
2.1.4.Etude du problème de contact 48
2.1.4.a. Conditions de contact 49
2.1.4.b. Régularisation de la loi de contact 50
2.1.4.c. Géométriede contact en 3D 51
2.2. Méthode Asymptotique N umérique 53
2.2.1.Méthodesde perturbation 53
2,2.2.Les approximants de Padé 56
2.2.3.Discrétisation par élémentsfinis 56
2.2.4.Procédure de continuation 58
2.3. Conclusions 59
-4t-
Chapitre 2 Formulation du problème de contact entre un coros rieide et une coque
Introduction
L'emboutissage, est un procédé de mise en forme des matériaux, généralementà froid. Il

consisteà déformer d'une façon permanentedes tôles minces sous I'effet des outils. Cette
définition se traduit numériquement,lors de la simulation, par une combinaisonde différentes
non-linéaritéstelles que la non-linéaritégéométrique(tôles minces),plasticité (déformation
permanente)et non-linêantê due au contact (outils). D'après l'étude bibliographique présenté
au chapitre 1, les lois élastoplastiques sont les lois de comportement qui décrivent
véritablement le comportementdes métaux au cours du procédé, et les élémentsfinis du type
coque semblentles plus efficacesdans la discrétisationdes tôles destinéespour la mise en
forme.
Dans ce chapitre, nous présenteronsune modélisation tridimensionnelle du contact

entre un corps rigide soumis à un déplacement aléatoire et une coque élastoplastique.
L'objectif est de proposerune simulationde I'emboutissageen utilisant des élémentsfinis du
brpe coque. Dans cette simulation, nous adopteronsune formulation de coque bien adaptée
aux phénomènesde grands déplacementset de grandes rotations. Au niveau de la loi de
comportement et sous I'hlpothèse que dans ce procédé les phénomènes de grandes
transformations (grands déplacements et grandes rotations) sont plus importants que les
phénomènes de grandes déformations, on va adopter une loi élasto-plastique en petites
déformations. Dans le cadre de la MAN, pour surrnonter la difficulté engendréepar la
discontinuitéde cetteloi, une procédurede régularisationest employée.
Le contact entre la tôle et les outils est assimilé à un contact entre un corps déformable
et des corps rigides. Ces conditions de contact constituentune très forte non-linéarité due à la
singularité de ces conditions. Comme pour la loi de comportement, pour appliquer une
technique de perturbation, une procédure de régularisation a été proposée dans le cadre des
Partantde ces travaux,nous allons
travaux de thèsesd'Elhage-Hussein[Elhage-Hussein/98].
faire une extension de la modélisation du contact à une échelle 3D et nous allons imposer des
déplacementsaux surfacesrigides pour décrire le mouvement des outils lors d'une opération
de mise en forme.
-42-
Chapitre 2 Fonnulation du problème de contact entre un corps rigide et une coque
L' objectif de ce chapitreest de présenterle problèmeà résoudre,d'une part la théorie

de coqueutilisée, la loi de comportement,les conditionsde contactet d'autre part le principe
variationnelqui lui est associé.
-43-
Chapitre 2 Fotmulation du nroblème de contact entre un corps rigide et une coque
2.1. Formulation du problème
Dans cette étude, nous considéronsune formulation tridimensionnelle du problème de

coque en présencedes conditions de contact unilatéral. La formulation adoptéeest baséesur
la formulation de coque proposée par Bûchter et Ramm I Bùchter et RammÆ2]. Une telle
formulation a été choisie pour deux raisons; premièrementla non-linéarité géométriquepeut
être exprimée dans un cadre quadratique ce qui est très conseillé dans la littérature des
M.A.N., en fait, elle permet une application facile des techniques de perturbation et
I'obtention d'un algorithmerapide et donc un gain en temps de calcul. Deuxièmement,cette
formulation utilise une loi de comportement tridimensionnelle sans condensation ce qui
permetune implémentationmoins difficile des lois de comportements
non-linéaires.
2.1.1.Description cinématiquede la coque
Dans cette formulation, la cinématique de base utilisée repose sur les théories
classiquesdesplaqueset coques.
Figure (2. I ) : Géométrie et cinëmatique de la coque
Ainsi un point quelconque de la géométrie de la coque (figure 2.1) est repéré, sur la
configuration initiale, par un vecteur x exprimé sous la forme suivante:
x10 t,e ',et) = r ( 0t,0' ) +0 3ar (e t,e' )
-44-
où r est un vecteur qui repère la surfacemoyenne, a3 est le vecteur directeur de la coque au
point considéréet (et,0',d'; désigneles coordonnéescurvilignesconvectives.En supposant
une variation linéaire du champ de déplacementu dans la direction de l'épaisseur, ce dernier

peut s'écrirecomme suit:
t
u 1 0 ' , 0 ' , e t ) =v ( 0 ' , 0 ' 7 +0 w ( 0 t , e ' )
où v et w représententrespectivementle déplacementd'un point de la surfacemoyenne et la

difference entre les vecteurs directeursde la coque entre la configuration de réferenceet celle
déformée.
2.1.2.Equation d'équilibre
Les inconnues du problème dans la formulation de coque adoptéesont le champ de

déplacement u, la déformation de Green Lagrange 1 et le secondtenseurde Piola Kirchhoff
S. En 3D l'équilibre du systèmes'exprimepar la relation suivante:

's,ayar-q(ôu):e (2.r)
f
où {(âu) est lè travail générépar les forces de contact que nous expliciterons ultérieurement.
Pour I'instant, on supposeque les forces de contact sont données,ce qui entraîne la linéarité
de la fonctionnelleP.(.).
Dans ce travail, nous utiliserons l'élément de coque proposé par Bûchter et Ramm
[Bi.ichter et RammÆ2] avec une variation linéaire du déplacementdans l'épaisseur comme

dans les formulations conventionnelles de coques. En revanche, il inclut une variable
supplémentairequi représenteune variation linéaire de la déformation ytz dans l'épaisseur.
Cette variable permet en effet, d'utiliser une loi de comportement tridimensionnelle sans
condensationmais aussi de prendre en compte les effets des grandesdéformations. Elle est
introduite au moyen du conceptE.A.S. de Simo et Rifaï [Simo et Rifai790].Le principe de ce
concept est d'introduire une déformation additionnelle I indépendantedu déplacementet qui
doit être orthogonaleau champ des contraintes:
1-Bu =l (2.2)
-45-
où Bu est la partie compatible du champ de déformation (Bu=lu)+y,,(u)). La
déformation additionnelle I doit vérifier la condition additionnelle suivante:
=o
l't,7 dv (2.3)
La présenteformulationn'est pas différentedu problèmeinitial : en fait, la solutionexacte,

correspondà
T=0 .
En élasticité,l'équation(2.3)découlede la stationnarité
de la fonctionnelle:
j J'Oo+D:D:@u+])dv (2.4)
où Dest la matricede comportement.
2.1.3.Loi de comportement
Au cours du procédé de l'emboutissage,les matériaux subissentdes déformations

permanentes.Dans cette section,nous tiendronscompte des déformationsplastiquesdans la
simulationdu procédé.
La loi de comportementen présencedes déformationsplastiquesdépendde I'histoire du
chargementd'où I'aspectincrémentalde la plasticité[Hi11/50][Halphenet Salençon/87].
Pour traiter la plasticitédansle cadrede la M.A.N., destravaux ont été menéspar Brunelot
[Brunelot/99] pour la résolution des problèmes viscoplastiquesen grandesdéformations.Un

autretravail est en cours, pour modéliser de vraies lois élastoplastiquesqui tiennent compte de
la déchargeélastiqueflmazaten/0ll.
Dans cette étude,la déchargeélastiquen'est pas prise en compte,nous adoptonsdonc

une loi de comportementplastique basée sur la théorie de déformation totale. Ce qui est très
convenable quand la non-linéarité physique est plus importante que l'effet du processus
irréversible et I'histoire du charsement.
La difficulté du traitementde ce tlpe de non-linéaritérésidedans I'aspectdiscontinudes lois

mises en jeu et la formulation non-quadratique.Pour obtenir une représentationanalyique et
pour que le problème soit bien posé, deux étapessont nécessaires.La première consisteà
régulariser la loi de comportement car cette dernière n'est pas analytique et n'est donc pas
développableen séries entières.La deuxième étape doit conduire, grâce à un changementde
-46-
Chapitre 2 Formulation du oroblème de contact entre un corps rigide et une coaue
variables, à réduire le degré de non-linéaité du problème en introduisant des variables

intermédiaires. Ce qui permet une écriture quadratique du problème et par conséquent
d'améliorer l'efficacité numérique de I'algorithme asymptotiqueen réduisant de manière
significative le temps de calcul. La taille mémoire nécessairepour le stockagedes variables
intermédiairesrestenégligeabledevantcelle du problèmeglobal fPotier-Ferryet all97].
Dans cette partie, nous avons choisi une loi de comportementbaséesur la relation de
qui s'écrit dansle cas 3D sousla forme suivante[ABAQUS/98-l] :
Ramberg-Osgood
-'
Ey=(r+ - (t- zù Pr+
v)so t' (2.s)
?l*]'
avec E, v) a, n,So,P,I et 6y sont des paramètresintrinsèquesdu matériauet
représententrespectivementle module de Young, le coefficient de Poisson, un coefficient

scalaire, I'exposant réel d'écrouissage,le tenseur déviateur des contraintes, la pression
hydrostatiqueéquivalente( P=jSI ), la matrice unité et la contraintelimite d'élasticité.
Comme le coefficient d'écrouissagen'est pas un entier, cette relation n'est pas analyique
pour des contraintes nulles, une procédure de régularisation a été introduite en modifiant
légèrementla contrainte équivalentede Mises ,S,c.
sr' jso,so+rfi,o, (2.6)
où r7. est un paramètrede régularisationpositif assezpetit. Remarquonsque pour 4. =0 nous
retrouvons la loi de comportementinitiale.

Une nouvelle variable est introduite sous la forme suivante :
-'
=+Â
* " &] =if *7,'l+
' L t o;-rs',s'
y (2.7)
Loyl I
Pour surmonter la difficulté due au terme en puissance(n est réel), dans le travail de Braikat
[Braikatig5] I'ordre de troncature des sériesa été limité à I'ordre 6. Zahrouni et al lZahrouni
et all98l ont évité cette limitation en utilisant une équation differentielle et en introduisant la
variablesuivante:
-47-
sii==l-So :So+(r7
ç' = fi' (2.8)
OY IOY
ce qui signifie que :

n-l
K =tûç
7 -
(2.e)
Cetterelation entre rc et Ç est une relation en puissancequi n'est pas quadratique.Toutefois
une fonction de type y=x'peut être écrite sous une forme differentielle équivalente
xdy=ty*, qui est quadratiquepar rapport aux variables x et y. Les deux variables rc et Ç
sontdéfiniesdonc par la relation suivante:
Çdrc=(n-r)rdÇ (2.10)
Pour obtenir la même pente au démarrageque celle de la loi sansrégularisation,il est facile
de montrerque le premier termede l'équation(2.5) doit êtremultiplié par (l+an) .
Donc la loi de comportementrégulariséese traduit par les équationssuivantes:
oE ( I + a ( 4 ^ ) ' ) y = ( l + v ) S ' - ( t - 2 v ) P I + r S '

^2
oÇ2 =è:s-= 3 So:So
+(r7fi' (2.rr)
Oy 2or'
oÇdrc=(n-l)rcdÇ
2.1.4.Etude du problème de contact
Dans la simulation numérique des procédésde mise en forme, la prise en compte des
conditions de contact revêt un caractèreessentielpour une bonne mise en æuvre et un suivi
effectif du procédé [Wang/8a] [Hacquin/96] [Chenotl9T]. Ces conditions de contact sont
évolutives au cours du processus,le matériau à déformer se rapproche des outils, prend le
contact,reste en contact (adhésion)puis il le perd (décollement).Au niveau modélisation,le
problèmede contactest gérépar un ensemblede conditionsunilatéralesfKikuchi et Oden/88]
[Raous et all88] [Jean/93]. Ces conditions constituent de nouvelles difficultés à cause de

I'aspectnon-linéairede ces conditions[Onâteet Saracibar/92].
Comme dans la modélisation de la non-linéarité matérielle, pour appliquer une

technique de perturbation au traitement du problème de contact, une procédure de
régularisation a été proposée par Elhage-Husseinet al lBlhage-Husseinet al 981. Le travail
-48-
Chapitre 2 Formulation du problème de contact entre un corps riside et une coque
avait pour but la résolution du problème de contact bidimensionnel sans frottement entre une
structureélastiqueet des solidesrigides.
Pour développerun modèle 3D utilisant des élémentsfinis de tlpe coque,nous avons
effectuer une extension de ces travaux à une échelle tridimensionnelle et nous avons tenue
comptedes mouvementsdes corpsrigides,représentantles outils, au coursdu processus.
2.1.4.a- Conditions de contact
Nous considéronsune structureélastoplastiquedont la loi de comportementsatisfait la

relation de Ramberg-Osgoodmentionnée précédemment.Nous supposonsque cette structure
occupe un domaine f) dont une partie de la frontière ôO" est susceptibled'entrer en contact
avecune surfacerigide.
Corpsrigide
Figure (2.2) : Contact : interaction surface rigide/surface dëformable
La zone effective de contact ainsi que les efforts de contact évoluent au cours du
processuset sont donc les inconnuesdu problèmeétudié.L'action localede corpsrigide sur la
structureest représentéepar une force de réaction notée R.
Un point P de la surface déformable entre en contact avec la surfacerigide lorsque la
distancequi les séparedevientnulle.
-49-
Chapitre 2 Formulation du problème de contact entre un corns rigide et une coque
Figure (2.3) : Loi analytique de contact
La vat'rable h est la distance ou l'écart entre les deux corps au point de contact
considéré, n est la normale extérieureà la surfacerigide passantpar ce point.
A tout point candidat au contact sont associéesles conditions classiquesdu contact

unilatéral suivantes:
o Le point de contact ne doit pas pénétrerla surfacerigide :
h>0 (2.r2)
o Le point de contact ne doit pas adhérerà la surfacede contact :
R.n)O (2.r3)
o D'une part lorsque le point est en contact, l'écart h est nul et une réaction
de contactest générée.D'autre part en I'absencede contact,h est non nul et
la force de contactest réduite àzéro;
(R.n).ft=0 (2.r4)
2.1.4.h.Régularisationde la loi de contact.
Les conditions de contact (2.12), (2.I3) et (2.14) ne sont pas analytiques.Or cette
condition est nécessairepour appliquer une technique de perturbation. Elhage-Husseinet al
[Elhage-Hussein et /98] ont proposé une régularisation par une relation hyperbolique, il ont
ainsi écrit:
-50-
Chapitre 2 Formulation du problème de contact entre un corps rigide et une coaue
*=dîEo (2.rs)
où â est l'écart initial entre les surfaceset 7 est un petit paramètrepositif suffisammentpetit
( I'indice c permet de I'identifier par rapport au paramètreutilisé pour la régularisation de la
loi de comportement(2.11) ). Ainsi pour une valeur assezpetite de r7", les conditions de
contact(2.12), (2.13) et (2.14) sontpresquevérifiées.
La surface de contact est constituée d'un ensemblede points caractériséschacun par

un couple de variables ô et h, cettedernière régularisationn'est donc pas uniforme sur toute
la toute la surface. Nous devons donc définir le paramètrede régularisation en chaque point
afin d'homogénéiser le comportement à I'interface. Si nous supposonsun écart donné,ho et
que nous lui associonsune pressionde contact &, alors le paramètrede régularisationest
défini en chaquepoint de contactpar I'expression suivante :
4,(point):#z (2.16)
Ainsi pour I'ensemble des points de contact, la régularisation est uniforme et caractérisée
par les deux paramètresho et Ro.Bien évidemment,lavaleur de hodoit être inférieure à la
valeur minimale des écartsinitiaux â de chaquepoint.
2.1.4,c. Géométrie du contact en 3D
Pour résoudrele problème,en plus de la relation (2.15), deux équationssupplémentaires

sont nécessaires.Ces deux équations sont données par la géométrie des surfaces rigides
considérées.Elhage-Hussein[Elhage-Hussein/98]a étudié dans un cadre bidimensionnelle
problème de contact pour différents types de surfaces, surface plane, circulaire et surface
arbitraire définie par une courbe de Bézier. Dans ce chapitre nous nous limiterons à l'étude de
contact entre une coque qui a un comportementélastoplastiqueet une surfacerigide sphérique
soumiseà un déplacementarbitraire. C'est ce déplacementqui va nous permettre de décrire le
mouvement des outils au cours du processusde I'emboutissage.Le traitement du problème de
contact en considérantdes corps rigides de formes complexesseratraité ultérieurement.
-51-
Corps rigide
C onfieuration in itiale
Confisuration finale
Figure (2.4) : Modélisation du problème de contact (rëférencesgéométriques):

cas d'un corps rigide soumis à un dëplacementarbitraire.
Dans la présente modélisation, la déformation dans la structure est due seulement au

déplacementdu corps rigide que nous avons choisi ici sous la forme de .Lw. Le vecteur
w représentele déplacementimposé sur I'outil, donc connu en module, senset direction. 2
est le paramètrede contrôle du problème. Un point P de la structureest défini par le vecteur
position X et les paramètresinitiaux de contact(ôno ) dansla configurationinitiale. Il occupe
la position P dans la configuration courante aprèsun déplacementu (figure 2.4). L'êcart de

contact h etla normale n peuvent donc être expriméspar ces deux équationsquadratiques:
(h+ p\'z=(X+u-C -2w, X+ u-C -2w ) (2.r7)

(h+p)n=X+u-C-Zw ( 2. 18)
où p, X et C représententrespectivementle rayon de la sphère , un vecteur décrivant la
position initiale du point étudié et un vecteur décrivant la position initiale du centre de la

sphère. Ces deux équations peuvent être substituéesdans l'équation (2.15) pour obtenir
I'expression de la force de contact R en fonction, seulement, du déplacement u et du
paramètrede charge2.
-52-
Dans ce travail, nous nous limiterons à une formulation de contact 3D sansfrottement.

La modélisation du problème de frottement dans le cadre de la M.A.N. n'a pas posé de
problèmesparticuliers dans les travaux de thèsede Brunelot I BruneloV99].
2.2. Méthode AsvmptotiqueNumérique
2.2.1.Techniquede perturbation
Dans cette section,la méthode asymptotiqueseraappliquéepour résoudrele problème

global correspondantà l'équation d'équilibre (2.1), la condition de compatibilitê (2.3),laloi
de comportement(2.11) et les conditionsde contact(2.15), (2.17) et (2.18) . Nous appliquons
donc une technique de perturbation qui consisteà développertoutesles variables du problème
I'ensemble
en sériesentières.Nous notonspar U=U(urZSrfr,nrR...) le vecteurreprésentant
des inconnues. Ce vecteur et le paramètre de charge sont développéssous forme de séries
entièrestronquéesà un ordrep au voisinaged'une solution initiale connue(U o,4) par rapport
à un paramètrede développement<<a >>:
o* a(À+a2(Jz+....,
fugl=u (2.re)
1"
I l(a1=1o* ah +a-k +.....
Le paramètre <<a > constitue une inconnue additionnelle que nous choisissonsde définir par
I'expressionsuivante:
a = ( u - u o , ur ) + Q " - ) , 0 ) , r (2.20)
Le choix de ce paramètreest basé sur la même stratégieque celle adoptéedans les méthodes
itératives. La forme de la dernière équation coffespond à une combinaison de déplacementet
de charge comme la longueur d'arc imposé dans le cas des algorithmesitératifs. Nous aurions
pu aussi identifier ce paramètre à un incrément de charge ou de déplacement [Cochelin et
all94l. En substituantles équations(2.I9) et (2.20) dans celles qui gouvernentle problème
global et en identifiant membre à membre par rapport à la puissancede <<a >, on obtient une
séquencede problèmes linéaires bien posés qui vont être résolus par une méthode de
discrétisationpar éléments finis jusqu'à un ordre de troncature de séries. Ces problèmes
-53-
linéaires admettent tous le même opérateur tangent et donc la même matrice de rigidité
tangente. Nous obtiendrons ainsi une représentationanalyique d'une partie de la branche
solution en n'effectuantqu'une seuledécompositionde la matrice de rigidité tangente,ce qui
est très important pour des problèmes fortement non-linéaires et surtout pour des problèmes
de grandetaille.
Le développementde l'équation d'équilibre et la condition de compatibilitéjusqu'à un ordre

de troncature p donne :
(2.2r)
[{'"rt r,,1uo,6u)+'s ,:2y,,(uo-,,6u1p"-
o:(y,(âu)+y,,(uo,ar;;+!'s [,^;*or=o
(2.22)
['s o:6Vdv=o
Le développementdes équationsconstitutives conduit au systèmesuivant :
(
lOrdre1: S,=D,:7,
I (z.zt)
p: So=D,:yr*Si'
lOrdre
où la matrice D représentele module tangentde la loi de comportementet qui a I'expression

suivante :
D,=k,s+crrr+c,sfsf
] ez4)
où I est le tenseur identité, 5 est un tenseur identité d'ordre quatre,Ct , Cz et G sont des
constantesdéfiniespar les expressionsuivantes:
n ._E(l+ary)
-'-
çyço
n^_(3 v+rc')Ct
wz- y1-2.v)
/-^- -KCr
"'-1*tra6oaKS{Sd
xJ$ilr"
z06ç6
-54-
Sf,'est le tenseur des contraintes résiduelles (résiduelle au sens de développement
asymptotique)à l'ordre p :
#rJc[r'a*a.rir'sÉ-')sd -o*o (2.25)

]-t].t E*,rr-,
C, et I sont desconstantesdéfinies par les expressionsuivantes
r.-- -K
-'-
1*travç'o
1rKSf,:sd
L'expression des contraintesest ensuite introduite dans l'équation d'équilibre pour

obtenir un problème avec le champ de déplacementcomme seul inconnue et par la suite une
formulation d'éléments finis en déplacement.
Après développement en séries entière et identification membre à membre par rapport à la

puissancede <<a >>,les conditions de contact(2.15),(2.17) et (2.18) donnentles expressions
suivantes:
A I'ordre1:
(X+uo-c-Zow , u,-2,w)
h,=
h o +P
tur-).w -hrn, (2.26)
nl =
h o +P
*,=ilîDo, - o+R o)
hL(r7"n
A I'ordrep :
hP
ho+p 2A hr+p
up-).,w-hpno $h4,-,
(2.27)
î, = -à-i,.p
o*,
-f;{, -.+R,-,)
",=o#n o "no**,>ifrt,"n,
Nous pouvons remarquer qu'après substitution des expressionsde l'écart hp et la normale
np dans celle de la réactionde contact Rc, cette dernièreva dépendreseulementdu champ
-55-
du déplacementet du paramètrede chargement)., ce qui va rendrefacile I'application de la

méthode des élémentsfinis.
2.2,2.Les approximants de Padé.
Dans la littérature des méthodes asymptotiquesnumériques,les représentationspar les séries

ont été amélioréesen introduisant les représentationspar les fractions rationnelles nommées :
approximantsde Padé.[Najah et aUl98l,[Elhage-Husseinet all00]:
I u,or=u o*o !-u,

-n_r I + o' l'l-.tJ2+...+d'-r:f-u
D,-t-z
"'
Dn r n
-
| (2.28)
D=n-'
f ^' - rt "o' r =' -^o ^' * * oD % -1,+az ^,*...+a'-' I I
n-r"t'.- Dn.r'"''""' Dn-r"n
|
on D@) sontdespolynômesde degré(i ) avecdescoefficientsréels@i)[i=r,n-r] :
D(a)=la sdr+ azdz+...+aidi (2.2e)
Il a été démontré dans [Najah et al/198], [Elhage-Husseinet al/001que les représentations

rationnelles augmentent considérablementle domaine de validité des séries et par suite la
taille despas de calcul ce qui permetune réductiondu tempsde calcul.
Dans ce qui suit , les résultatsprésentéscorrespondrontà un algorithmede la M.A.N. avecles
approximations de Padé. En fait dans ce travail, cette technique a été implémentée et a
montré, via destests,une améliorationdes domainesde validité des séries.Généralement,les
approximants de Padé permettent une réduction de I'ordre de 50Yo du nombre des pas
pour trouver la solution.
nécessaires
2.2.3.Discrétisation par élémentsfinis
Nous utilisons une discrétisation usuelle d'un élément quadrilatéral isoparamétriqueà

huit næuds. La modélisation de la géométrie et du champ de déplacement s'expriment
respectivementpar:
-2.
a=)N*rt + 0t Noa\
k=l
4
u=)Ntvt+0tNkwk
k=l
-56
Chapite 2 Formulation du problème de contact entre un coros rigide et une coaue
où Nk désigne les fonctions de forme, 1o, a\, vk et wt sont des variables nodales
interpolées via les fonctions de formes classiquesdes élémentsquadrilatéraux à huit næuds,

k indique le numéro du næud et z désignele nombre des næudsde l'élément qui vaut 8 dans
l'élément fini considérédanscetteétude.
Après I'application des techniques de perturbation et la discrétisation par éléments

finis, le problèmeà I'ordrep peut s'écrire sousla forme suivante:
lr,\uo\=V;'
l.t^;l (2.30)
où [r(,] est la matrice de rigidité tangentesans contact,elle est calculéeau voisinagede la
discrétiséà l'ordrep, V;'} est le

solutioninitiale (Uo,,a), {Uo} est le vecteurdéplacement
second membre classique dans la méthode asymptotique numérique. Le vecteur {Uo } est
desnæudscandidats
en deuxtermes' 1ur)={ui ,Uïl,où tÛ} estle déplacement
décomposé
|
au contactet lue" "jt le déplacement Le vecteurt";Fh;,Ol
du restedesnæudsde la structure.
les forcesde contactet nouspouvonsI'exprimersousla formesuivante:
représente
îl*uft"
t";Ft""ttu l.b;'l (2.3r)
[K"] est une matricesymétriquequ'on appellematricede rigidité de contact.$"] .rt ort

vecteur force qui dépend surtout du vecteur déplacementimposé ,n . {F;" } est ,rn vecteur qui
dépendnon-linéairementdestermes de contact aux ordresprécédents.
Le problème est discrétiséen utilisant l' élément fini coque mentionnéprécédemment.Pour

assurerune bonne description des phénomènesde flexion, 5 points d'intégration numérique
de Gausssont utilisés suivantl'épaisseur.
-57 -
2.2.4.Procédure de continuation
Comme les séries ont un domaine de validité limité, une procédurede continuation a
été proposéepar Cochelin [CocheIin/94]pour décrire toute la solution.Dans chaquepas la
valeur maximale (a^) du paramètre <<a >>est définie automatiquement.Un simple critère
consiste à exiger que la solution entre deux ordres consécutifs reste inférieure à une valeur
a'.
I
I ll''ll)F
o'=["iii (232)
)
Dans ce travail, nous utiliserons une procédure de continuation basée sur les
approximants de Padé qui a été proposée récemment par Elhage-Hussein et al [Elhage-
Hussein et all00]. Dans ce cas, le domaine de validité (a^p) de la solution est défini par
I' expressionsuivante:
llp,@*)-P,-,(aQll
6'=-1iF,(o;)_u;lj-- (2.33)
oir la valeur maximale anp est cherchéedans I'intewalle la^, premier pôtel et déterminée
par une technique de dichotomie. Le premier pôle est la plus petite racine réelle du
dénominatenrDnt de l'équation (2.29).
-58-
Chapitre 2 Fonnulation du problème de contact entre un corps rigide et une coque
2.3. Conclusions
Dans ce chapitre, nous avons présentéla formulation 3D d'un problème de contact entre une
coque élastoplastiqueet une sphèrerigide soumiseà un déplacementarbitraire. Les conditions
de contact évolutives sont transforméesen une relation analytique fonction de l'écart entre la
structure et la surface rigide. Il s'agit donc d'une loi de contact régularisée.Ceci nous
permettra de considérer un ensemble des points de la structure, susceptibles d'entrer en
contactavecle corps rigide. La résolutiondu problèmeva donnerl'évolution de cespoints de
contact,l'évolution de l'état de la structuredéforméesousI'effet du déplacementimposé à la
surfacerigide.
Après avoir formuler le problème dans un cadre quadratique, nous avons appliqué une
technique de perturbation. Ceci nous a permis de transformer, comme lors de tout traitement
d'un problèmemécaniqueavec la M.A.N, les problèmesnon-linéairesde départen une suite
récurrentede problèmes linéaires bien poséset admettanttous le même opérateurtangent.Ce
dernier contient des termesprovenant de la non-linéarité matérielle et des termesprovenantdu
contact. La méthode des éléments finis a été appliquée pour effectuer une discrétisation
spatialedes problèmeslinéairesobtenus.Ces problèmeslinéairesdiscrétisés,admettenttous
la même matricede rigidité, ce qui constitueun despoint forts de la méthodeutilisée(MAN).
Nous avons obtenu, donc un algorithme qui nous permettra de simuler numériquement
l'emboutissageen 3D des coquesminces où la non-linéaritégéométriqueest coupléeavec la
non-linéaritématérielleet la non-linéaritédesconditionsévolutivesde contact.
-59-
Chapitre 3 Application à la simulation de I'emboutissasehémisphëriaue
Chapitre 3
Application à la simulation de I'emboutissagehémisphérique
Introduction 61
3.1.Simulation de I'emboutissagehémisphérique 6l
3.1.1.Présentationdu problème 6l
3.1.2.Traitement desproblèmesnon-linéairesdans Abaqus 64
3.1.3.Comparaisondes résultats 67
3.2.Etude de I'ordre optimal de troncature desséries 69
3.3. Etude de l'influence des paramètres de régularisation 7l

3.3.1.Influence de la régularisation de la loi de comportement 7l
3.3.2.Influencede la régularisationde la loi de contact 72
3.4. Efficacité de I'algorithme 75
3.5.Etude de I'amincissementde la tôle 77
3.6.Etude de la qualité de la solution 79
3.7.Conclusions 81
-60
Chaoitre 3 Application à la simulation de I'emboutissagehémisphërique
Introduction
Dans ce chapitre, nous appliquons la formulation présentéedans le chapitre précédent

(chapitre 2) à la simulation de I'emboutissagehémisphérique.Ce procédé est considéré
comme un test standarddansI'industrie de la mise en forme desmatériaux.C'est un test qui a
été étudié expérimentalementpar [Gosh et Hecker/75] et traité numériquement pour la
première fois par V/ang et Budiansky fWang et BudianskylTS] dans un modèle basé sur une
formulation d'élémentsfinis de type membrane.Ce test est souventutilisé pour valider les
modèleset les simulationsnumériquesde l'emboutissage[Bellet et aV90][ABAQUS/98-2].A
travers la simulation de ce procédé, nous allons valider I'algorithme asymptotique présenté
dans le chapitreprécédent(chapitre2).
Les résultats obtenus seront comparés avec ceux du code industriel Abaqus, que nous
utilisons ici commeréférence.
3.1. Simulation de I'emboutissagehémisphériqued'une tôle circulaire
3.1.1.Présentationdu problème
L'emboutissage hémisphérique consiste à déformer une structure, généralementde

forme circulaire, par contact avecun outil (poinçon) de forme hémisphérique.
Sur la figure (3.1), on montre la géométrieet les propriétésmécaniquesde la structure
étudiée.La tôle est de I mm d'épaisseur,de forme circulaire et supposéeencastréesur son
contour. En raison de la symétrie de la structureconsidérée,seulementun quart de la tôle est
étudié.Le matériaude la tôle considéréeest I'acier désoxydéà I'aluminium (aluminium-killed
steel ), supposé satisfaire la relation contrainte déformation élastoplastiquede Ramberg-
Osgood(équation2..5 duchapitre2).
-61
Chapitre 3 Application à la simulation de I'emboutissase hémisphériaue
Geométrie: Propriétés
Mecanique:
(w
R:59.18rnrn E:206.8Gpa
950.8 rrrn \80.3
ælrnrn ov=170MPa
r4.75
Y,J
'emboutissage
Figure(3 . I ) : Simulation de I hémisphérique: géométrie et propriétés
mécaniques
La tôle est déforméepar contact avec un outil sphériquesoumis à un déplacement

vertical qu'on a noté w dansle paragraphe2.1.4.c.du chapitre2.
La structure est discrétiséeen utilisant l'élément fini décrit dans le chapitre précédent
[Btichter et RammÆ2]. C'est un élément quadrilatéral à 8 næuds ( 6 degrés de liberté par

næud ) qui est utilisé ici avec 5 couchesde points de Gaussd'intégration numérique dans
l'épaisseuret 4 points d'intégrationpar couche.
Dans cette étude,on supposeque tous les næudsde la structuresont susceptiblesd'entrer en
contactavec la surfacerigide exceptésles næudsdu contour qui sont encastrés.Par suite, le
nombre des points de contact considérédans les calculs et les discussionsconespond au
nombre de degrés de libertés du problème divisé par 6 ( nombre de degrés de libertés par
næud). Donc, un problèmede 6000 d.d.l., contientenviron 1000points de contact.
-62-
Chonitrc 7 Annlirnlion à ln çimulntinn dc l'emhoutissnçp hémisnhérioue
Figure(j.2) : Déformé de la tôle après un déplacementvertical

du poinçon de 27.8 mm
La figure (3.2 ) montre la déformée de la tôle qui coffespond à un déplacementvertical de

I'outil de 27.8 mm. Dans toute l'étude qui suit, le déplacementverticalede I'outil est toujours
limité à cette valeur maximum.
Au chapitre précédent,nous avons vu que les conditions non analytiquesde la loi contact ont
été régulariséespar une loi hlperbolique fonction de l'écart entre le point de contact et la
surfacerigide (h) et le paramètrede régularisation4. Plus le paramètrede régularisationest
petit ( Rd est petit ), plus la loi de contact régularisée se rapproche du contact parfait. La
figure (3.3) illustre ceci.
-63-
Chapitre 3 Aoolication à la simulation de I'emboutissase hémisphériaue
Figure(3.3)VariationdeI'intensitédelaréactiondecontactenunnæudenfonctiondel'écarth ;
influence du paramètre de régularisation
La figure (3.3) représentel'évolution de la réaction de contact, en un point de la surface

considéréeen contact lorsque l'écart h entre ce point et la surfacerigide tend vers zéro. Pour
Rd=0.5, le paramètrede régularisationcorrespondantest très petit et la loi de régularisation
( À=r7
correspondante ôll
,,h est très voisine de la loi discontinue( Rft=O).
3.1,2.Traitement des problèmesnon-linéairesdans Abaqus
Fondé en 1978, Abaqus est un logiciel de simulation pour I'analyse de nombreux

problèmesd'ingénierie.C'est un programmegénérald'élémentsfinis qui offre la possibilité
de traiter et de créer différents modèles en combinant differentes procédures,différents types
d'éléments et différents matériaux. Pour traiter les non-linéarités de tlpe matérielles, Abaqus
permet de modéliser la plasticité des métaux avec des modèles de Drucker-Prager,la plasticité
avec état critique, plasticité des moussesécrasables(absorption d'énergie), plasticité des
métaux poreux, modèles de roches fissurées, influence de la vitesse de déformation,
viscoplasticitépar fluage etlou par gonflement,critère de plasticité de Hill et le modèle de
Romberg-Osgooden utilisant la théoriede la plasticitéen déformationtotale.
-64-
Chapitre 3 Aoplication à la simulation de I'emboutissaee hémisphériaue
Numériquement,le traitement des problèmesnon-linéairesdansAbaqus est basésur la

méthodede Newton-Raphsonavec une longueur de pas adaptative.
Pour traiter le problème de contact, Abaqus utilise en général la méthode des
multiplicateursde Lagrange(5.4.d.du chapitrel).
Dans ce code,deux lois de contactsont disponibles:
la loi rigide ou dure: aveccette loi, deux surfacessont en contactsi la distancequi les sépare
est nulle. Une pression de contact est alors transmise à I'interface. Lorsque cette pression
atteint une valeur nulle, les surfacesse séparent.
p=0 si ft>O
p=kh aveck infini, si ft=O
pressionde contact
Toute valeur de pression

.-- possiblequand il y a contact
pas de pressionquand il n'y a

pas contact
Ecartement
Figure (3.4) Loi de contact non-régularisëe
Le contactde tlpe dur est la loi utilisée par défaut dans tous les problèmesde contactdans
Abaqus.
La loi régularisée: Dans ce modèle, la pressionest définie en fonction de l'écart (ft ) par une
relation exponentielle.L'utilisateur doit entrer une distanceinitiale de contact ( c ) et une
valeur de pression(po) qui correspondà la pressionen ( h=0). En effet, le contactentre les
deux surfacesest établi si (fr ) atteint la valeur (c).
-65-
Chapitre 3 Apolication à la simulation de l'emboutissase hémisohériaue
P=0 pour h<-c
o--ff!.t)(n'*."-l)] pour -cch<6c
o--fiha-r*çL-a76"';1f
pow 6c6cch
Pressionde contact
Loi exponentielle
Figure (3.5) La loi régularisée : courbe pression enfonction de h
L' algorith m e m aster-sI ave.
Lorsque I'approche'Contact Pair' est considérée,Abaqus utilise le concept 'master-

slave' pour définir I'interface de contact. Le principe de l'approche 'Contact Pair' est
d'associerdeux surfacepotentiellementen contact,de telle sorteque I'une des deux surfaces
ne puisse pas pénétrerI'autre. Une des surfacesest désignéesurface'esclave'et I'autre la
surface 'maître'. Les næuds de la surface esclavesont contraints de ne pas pénétrer la surface
'maître'.Les næudsde la surface'maître'peuvent,en principe,pénétrerla surface'esclave'.En
général, la surface 'maître' est choisie cofirme étant la surface du matériau le plus rigide si les
matériaux sont differents ou la surface du corps qui possède le maillage le plus fin si les
matériaux sont les mêmes. En présencede corps rigides, la surfacerigide est toujours la
surface'maître'.
-66-
Chapitre 3 Application à Ia simulation de I'emboutissase hémisphérique
Les options 'Small Sliding' et 'Finite Sliding'
Avec I'approche 'contact pair', deux formulations sont disponibles pour modéliser
I'interactionau niveau de I'interfacede contact.La premièreest la formulation'small-sliding'
dans laquelle les deux surfacespeuvent faire des rotations et de faibles glissementsI'une par
rapport à I'autre.Ceci implique que les næudsde la surfaceesclaveinteragissenttoujours avec
la même zone de la surface'maître'.La deuxièmeest la formulation'finite-sliding'qui est la
plus générale et qui est utilisée par défaut. Dans cette dernière, de grands glissements,de
grandes rotations, et de grandes déformations des surfaces sont autorisés. Pour éviter les
problèmes de convergence,Abaqus lisse automatiquementles surfaces maîtressesdans un
contact de tlpe 'finite-sliding'. De toute évidence la formulation 'finite-sliding' est plus
coûteuseen tempsde calcul surtouten 3D.
Représentation des surfaces de contact
Dans l'étude de contactentre des corps déformables,il est évidentque les deux corps
doivent être maillées.Pour représenterles surfacesde contact dans le cas des corps rigides,
trois approchessont possiblesdansAbaqus:
- Définir les surfacesanalytiquement,en 2D ou 3D
- Mailler la surfaced'un corpsen utilisant desélémentsrigides
- Ecrire un sousprogrammeutilisateur (RSURF) pour définir la surface
Dans ce travail, nous avons choisi d'utiliser la loi de contact 'rigide' non-régularisée.Nous
utilisons aussi I'approche'Contact Pair' qui est adoptéepar Abaqus pour traiter la plupart des
problèmesde contact [Abaqus/98-l]. Dans cette approche,on décrit uniquementles surfaces
en contact par une paire de surfaces,sansdéfinir explicitement les éléments de contact. Ces
derniers sont directement gérés par cette approche.Le corps rigide utilisé dans le problème
étudié dans ce chapitre est une sphèrequ'on représenteanalyiquement dans le calcul Abaqus.
3.1.3.Comparaisondes résultats
La figure (3.6) montre une comparaisonde la variation descontrainteséquivalentesde

Mises en fonction du paramètrede contrôle )" qui contrôle le déplacementdu poinçon. Les
courbescorrespondentau point de Gaussle plus prochedu centrede la tôle.
-67 -
Chapitre 3 Application à la simulation de I'emboutissase hémisohériaue
Figure (3.6) : Evolution de la contrainte équivalentede Mises enfonction du paramètre

de charge ), Comparaison avecAbaqus
On remarque la bonne précision de I'algorithme asymptotique pour estimer les

contraintesdansla tôle. Les courbessont décaléesvers la droite à causede l'écart initial entre
la tôle et le poinçon, qui vaut 2 mm danstous les calculs effectuésdanscette étude.
'o ?s 30
o,ll*.*,0".'l'-,*,
Figure (3.7) : Evolution de laforce du poinçon enfonction de son déplacement:
comparaison avec le code industriel Abaqus
La figure (3.7) représentel'évolution de la force du poinçon en fonction de son déplacement

au cours du processus.Cette force correspondà la force totale de contact qui est équivalenteà
-68-
Chapitre 3 Apolication à la simulation de I'emboutissaee hémisphériaue
la somme des forces de contact en chaque point candidat au contact. Dans les procédés de
mise en forme, cette force correspondà la force exercéepar les outils. Cette demière est aussi
appeléeforce d' emboutissage.
On peut remarquer, à partir de ces deux dernières figures, la bonne estimation des
contraintesdansla tôle ainsi que la très bonneprévision de la force de I'emboutissage.
3.2.Etude de I'ordre optimal de troncature des séries
Nous avons vu dans la revue bibliographique et le chapitre 2 que I'algorithme asymptotique
fait intervenir un secondmembre {r''} dans les équationsqui gouvernentle problème.Le
tempsde calcul de ce terme dépendsurtout de I'ordre des séries.A partir d'un certain ordre,
ce temps de calcul devient élevé. Dans cette section, on va étudier le temps CPU afin de
déterminerI'ordre optimal de troncaturedes séries.Aucune comparaisondu temps de calcul
avec les méthodesitératives ne seraprésentéedans un premier temps. Cette comparaisonfera
I'objet d'une sectionultérieure.L'étude est faite pour un algorithmeasymptotiquebasésur les
représentationspolynomiales et un algorithme basé sur les représentations rationnelles
(Approximants de Padé) et ce pour différentes tailles du problème. Pour chaque taille du
problème et chaque ordre considéré,on présente,le nombre de pas nécessairespour achever
le calcul, le rapport entre le temps nécessairepour l'évaluation du secondmembre {F'} ., f.
pour l'évaluationet la décomposition

tempsft nécessaire de la matricede rigidité tangente
fX,l; "t finalementle tempstotal relatif correspondant

au rapportentrele tempstotal CPU et
le temps/t .
Ordre Nbre de pas (Fnl)/t(K0 tempstotal relatif

8 187 0,53 286,1l
l0 145 0,66 240.7
t2 123 0,86 228.78
t3 ll6 0,96 227.36
t4 ll0 1,08 228.8
r5 106 1,19 232.14
20 99 1,44 241,56
Tableau(3.1) : Ordre optimal,sëries:600d.d.l
-69-
Chapitre 3 Aoolication à la simulation de I'emboutissaeehémisohérique
Ordre Nbre de pas (Fnl)/t(Kt) tempstotal relatif

8 67 o5, 2s3.84
0 29 0,69 218,01
2 08 0.84 r98.72
J 01 0,94 195,94
4 97 1,13 206,61
5 93 t,24 208,32
20 80 r,92 233,6
Tableau (3 2) ' Ordre optimal, séries :1537 d.d.l
Ordre Nbre de pas t(Fnl)/t(Kt) temps total relatif

8 155 0-57 238,l 5
t0 ll9 0.67 98,73
t2 99 0.93 90,92
13 94 02 89,51
t4 90 t4 92,24
15 86 24 92.64
20 77 76 212.52
Ordre Nbre de pas (Fnl)/t(Kù Temps total relatif

8 ll0 0.49 163,9
l0 85 0.66 l4l.l
t2 IJ 0.85 135,05
13 68 0,95 132,6
t4 64 1.08 r33.12
l5 6l l.l9 r 33.59
20 56 1.43 136.08
Tableau (3 4) Ordre optimal, Padë :600 d.d.l
Ordre Nbre de pas (Fnl)/t(KO Temps total relatif

8 94 0.52 142,88
l0 75 0.7 t27.5
t2 63 0,84 r15.92
l3 59 0,94 t14.46
t4 )l r,04 116.28
l5 54 r,23 120,42
20 44 l.88 126.72
-70-
Chapitre 3 Anplication à la simulation de I'emboutissaeehëmisphëriaue
Ordre Nbre de pas t(Fnl)/t(Kt) Temns total rélatif

8 86 0,54 32,85
l0 71 0,66 t7,86
t2 60 0,93 15.57
t3 55 l.0l 10.56
l4 53 1.14 13.31
l5 5l 1.23 l 13,8
20 43 l.83 121.69
Tableau(3.6) : Ordre optimal,Padé:5766d.d.l
On peut ainsi observerla diminution du nombre des pas au fur et à mesureque I'ordre
de troncaturedes sériesaugmente.Un passagede l'ordre 8 à I'ordre 20 divise presquepar
deux le nombre de décompositions des matrices nécessairespour obtenir la solution.
Toutefois,comme I'ordre a augmenté,le tempsde calcul des vecteursseconds-"-br., {f"}

a aussiaugmenté:ce passagede I'ordre 8 à I'ordre 20 entraîneaussiune multiplication par
tF't
trois du rapport - . Comme le montrent ces tableaux, I'ordre optimal de troncature des
/K,
sériesest de 13. Ce résultatest en bon accord avec quelquerésultatsde la littératurede La
M.A.N. où l'ordre optimal varie souvent entre 10 et 15 et le rapport # utoirine I

/l(,
généralement.
3.3. Etude de I'influence des paramètres de régularisation
Dans la formulation du problème (chapitre 2), on a introduit deux paramètres de

régularisation.Il s'agit d'un paramètrede régularisationde la loi de comportementet d'un
paramètre de régularisationde la loi de contact. Dans cette section, nous allons étudier
I'influence de ces paramètressur le temps de calcul ainsi que sur la qualité de la solution
obtenuepar I'algorithme.
3.3.1.Influence de la régularisationde la loi de comportement
Pour étudier I'influence du paramètrede régularisatiofl 4^ de la loi de comportement
(équation2.9 du chapitre2), on a effectuéquatrecalculscorrespondants

à r7^=0.001, 0.01,
0.1 et ry^=|. Nous remarqueronsque pour les trois premièresvaleurs, I'algorithme donne la
-71
Chaoitre 3 Aoplication à la simulation de I'emboutissase hémisphérique
même qualité de solution et ce avec le même nombre de décompositionsde matrices de

rigidité tangente.Une valeur T-=1, induit une modification inacceptablede la réponse
globalefigure (3.8). Ceci paraît être tout à fait normal car dansce casla loi de comportement
a été énormémentmodifiée. Par la suite et pour tout le restedes calculs,nous prendronsune
valeurde 17*=0.1.
0.2 0.4 0.6 0.8

À
Figure (3.8) : Etude de l'influence du paramètre de régularisation

de la loi de comportemeilqn sur la qualité de la solution.
3.3.2.Influence de la régularisationde la loi de contact
Si le paramètrede régularisationde la loi de comportementpeut être choisi d'une

façon plus ou moins arbitraire dans un large intervalle, la régularisation des conditions de
contactdoit être définie avec beaucoupplus d'attention.Dans la loi de contact,pour défrnir
une régularisationforte ou une douce, on doit faire varier le couple (Rd,hd) de la relation
(2.16). Pour la même densité du maillage ( discrétisationen 600 d.d.l. ), on compare les
résultatspour Ra=5 et Ra=0.5 avecceux obtenuspar le codeindustrielAbaqus( rRd=O).
-72-
Chapitre 3 Application à la simulation de l'emboutissaqe hémisphérique
On remarque que le paramètre de régularisation de la loi de contact a une influence

considérablesur la qualitéde la solution.
e 600
'i soo
- - N,r.^.N.(Rd'0.5)
-0 4oo --- --- i\1.i\.N.{Rd,,,,,5}
Abaqus(Rd=0t
3 300
6 zoo
0.4
T
Figure (3.9) : Etude de l'influence du paramètre

de régularisation de la loi de contact Qcsur la qualité de la solution.
Comme le montre la figure (3.9) plus petite est la valeur de Ra, meilleureest la qualité de la
solution. Ceci n'est pas le cas pour la force du poinçon correspondante
à la force totale de
contact qui n'est pas très affectée par cette régularisation. En effet, au contraire des forces
locales (contraintespar exemple), les forces globales sont directementliées à l'équilibre
global qui est satisfait même en choisissant des valeurs très grandes du paramètre de
régularisationdes conditionsde contact(figure 3.10).
-73-
Chaoitre 3 Application à la simulation de I'emboutissase hémisphériaue
------ Ab3qui(Rd""0)
- - 1v1..{.N.(Rd-.0.5)
- Nl.^.N.(Rd''5)
l0 15 20
Déplacenentdu poinçon(m)
Figure (3.10) ' Comparaisonde laforce du poinçon:force d'emboutissage,

influence de la régularisation de la loi de contact
Maintenant, nous allons étudier I'influence de ce paramètre de régularisation sur le

tempsde calcul et par la suite sur la rapiditéde I'algorithme.
MAN.
Abaqus
RH RG{,5
185 39 68
Tableau (3.7) Influence de la régularisation de la loi de contact (hd:0.05), comparaison avec Abaqus :
Nombre de décompositionsdes matrices nëcessairespour atteindre un déplacementde 27.8 mm
pour un problème de 600 ddl
Dans le tableau3.7, on montre I'influence de la régularisationde la loi de contactsur

la rapidité de I'algorithme, comparé au code industriel Abaqus. On remarqueque le nombre
des pas de la M.A.N., correspondantau nombre de décompositionsde matrice,augmentepour
des petites valeurs du paramètrede régularisationdes conditions de contact.En effet, quand le
paramètre de régularisation est faible, la non-linéarité du problème devient plus sévèreet par
conséquenceplusieurspas sont nécessaires
pour acheverle calcul. On retiendrapar la suite
une valeur de Ra=0.5 (hd:0.05) pour assurerune précision suffisantedans I'estimationdes

variableslocales.
-74-
Chapitre 3 Application à la simulation de I'emboutissaee hémisphërique
3.4. Efficacité de I'alsorithme
Le problème présenté précédemment a été étudié pour trois tailles differentes du

problème,correspondantes
respectivementà 600, 1537et 5766 degrêsde liberté.
Dans le tableau (3.8), on compare le nombre des décompositionsdes matrices de rigidité
pour atteindreun déplacementdu poinçon de 27.8 mm. Rappelonsqu'ici
tangentenécessaires
la comparaisonest faite en terme de 'constructionset décompositions'de la matrice de
rigidité tangente.En effet, dans les méthodesitératives (méthodesutilisé par Abaqus, pris ici
comme réference), il y' a souvent plusieurs itérations à chaque pas et c'est le nombre de
traitement (constructions* décompositions) qui est déterminantpour le temps de calcul.
Comparé à Abaqus, l'algorithme asymptotique nécessite très peu d'opérations

'constructions et décompositions'. On remarque bien que le raffinement du maillage rend
difficile la convergencede la méthode itérative (Abaqus) à I'opposé de I'algorithme

'constructions-
asymptotiquequi préserveson efficacité. Abaqus nécessitesix fois plus de
décompositions'par rapport à I'algorithme asymptotiquepour un problèmede 5766 d.d.l.
d.d.l M.A.N. Abaqus

600 68 185
r537 59 t97
5766 55 326
Figure (3.8) Nombre de décompositionsde matrices de rigidité tangente :

comparaison de la M.A.N. avec Abaqus.
Notons au passageque pour les problèmes à grand nombre de degrésde liberté, c'est la
décomposition de la matrice de rigidité qui consommele plus de temps CPU, et que c'est sur
ce point que les efforts des chercheurssont particulièrementdéployés.
Notons que pour des problèmes de grande taille, le nombre de points de contact devient très
important. Ceci induit un mauvais conditionnement de la matrice de rigidité tangente qui
génèreune lenteurde la convergence.Par ailleurs,la taille de cettematrice dans les méthodes
itératives, qui utilisent les multiplicateurs de Lagrange(cas de Abaqus) dans le traitement du
-75-
Chapitre 3 Application à la simulation de I'emboutissaeehémisphërique
problème de contact, est plus grande que celle de I'algorithme asymptotique(chapitre 1,

5.4.b).Si on note par Nc le nombredespoints de contact,dansla méthodedesmultiplicateurs
de Lagrange,la matrice à inversera une taille égaleà (nddl+Nc)* (nddl+Nc), par contre dans
I'algorithme asymptotiquela taille de cettematriceest seulementde (nddl)*(nddl).Dans notre
algorithme,pour prendreen compteles conditionsde contact,on ne modifie pas la taille de la
matrice de rigidité tangente. Ce résultat est très important surtout pour le traitement des
problèmes de simulation de mise en forme où plusieurs non-linéaritéssont couplées,les
structures étudiées sont généralementde grande taille et un maillage très fin est souvent
pour assurerla convergenceet obtenirune solutionassezprécise.
nécessaire
On note que ces résultatssont obtenusen utilisant le même couplede régularisationde

la loi de contact(Rd,hd ). Si on analysele tableau(3.8), on remarqueune légèreréduction du
nombre des décompositionsdes matrices quand le nombre des points de contact augmente.
Ceci est dû au fait qu'aprèsdiscrétisation,chaquepoint de la surfacede contactest soumisà
une force R". Après raffinement du maillage, la force globale de contact reste constante
tandis que la force nodale de contact est réduite. Ce phénomène,comme le montre la figure
(3.9), est équivalent à une augmentationde la valeur du paramètrede régularisation de la loi
de contact ce qui adoucit sa non-linéarité.
â 600
.9 5m
c --
€ 400 ivtr\N 60ûtdl
.: l\t^N 1536ddl
-----M.4.N.5766d.d.|
'E
6 zoo
;t'
il
tl
o.2
Figure(3.1I) : Influence du maillage sur la qualité de la solution
-76-
Chapitre 3 Application à la simulation de I'emboutissagehémisphérioue
3.5. Etude de I'amincissementde la tôle
Parmi les phénomènesles plus importants qui se déroulent lors d'un processusde la
mise en forme des tôles minces, on rencontre le phénomènede I'amincissementdu flan, qui
peut engendrer des fissurations voire des ruptures. L'amincissement ou la variation de
l'épaisseur est parmi les facteurs qui doivent être bien estimés lors des simulations du
procédé.Il est I'un des facteursqui permettentd'évaluerl'emboutissabilitéd'une pièce dansla
phasedu bureaud'études.Le flan étudiédanscettesectiona une épaisseurinitiale de 1 mm.
Dans cette section, nous testonsnotre algorithme dans la prédiction de la variation de

l'épaisseuren le comparantavec le code Abaqus.
Etat initiale Etat Déformé
figure (3.12) Calcul de la variation de l'épaisseur
D'après la cinématiquede l'élément fini utilisé (figure 3.12),la valeur de l'épaisseure de la

tôle en un point donné se déduit directementdu vecteur directeur de la coque en ce point:
ll= ll
"=lloll
Ce vecteur directeur est calculé en fonction de I'orientation initiale de la coque et de sa
variation :
a3=a3+w
si on attribut les composantessuivantesà ces vecteurs ,u,1Ï,),^,,,,,, /)

vl
w)
1e,,e,,e,y
77-
Chapitre 3 Apolication à la simulation de I'emboutissase hémisohériaue
L'épaisseurde la tôle vaut donc :
+ull=,,!
e=lla, Q+tol +Q+,of+(r+.ol
Pour étudier la variation de l'épaisseur, on a pris un exemple correspondant à une

discrétisation en 600 d.d.l. du problème présenté précédemment.Le déplacementvertical
imposé au poinçon vaut toujours 27.8 mm. On compare à la fin du calcul la valeur de
l'épaisseurau niveau des næuds situéssur une sectionradiale.Au niveau du code industriel
Abaqus on a utilisé plusieurs tlpes d'éléments dont le calcul permet une évaluation de
l'épaisseur et sa récupérationdans les fichiers de sortie. La raison pour laquelle plusieurs
tlpes d'éléments sont pris est d'éviter I'effet de l'élément choisi dans le calcul. On a
remarqué (figure 3.13) une différence de I'ordre de 15 Yo entre les valeurs donnéespar
I'algorithme asymptotiqueet celle du codeAbaqus.
1.1
^ 1
(l)
0.9
<o
0.8
(.)
0.7 " M.A.N. (Ramrn)
(.)
q
- -- Abaqus(CAXI)
ô
0.6 Abaqus(53R)
rYl ----Abaqus(SaR)
0.5 " " Abaqus(C]AXS)
0.4
s t$uméro 20 2s
desnà3u4.
Figure (3. I 3) Etude de la variation de l'épaisseur de la tôle:

Comparaison avec le code Abaqus
En première hypothèse,nous avons supposéque l'algorithme de contact peut être à

I'origine de différence entre les résultats. Pour cela on a effectué un petit test de traction
axiale qui ne fait pas intervenir les condition de contact.La plaque étudiéeest discrétiséeen
utilisant deux élémentsde type Ramm. En premier lieu, nous avonsmené le calcul jusqu'à un
palier de déformation de 2%opour atteindre enfin une valeur de 20oÂ.Nous avons constatéque
-78-
Chaoitre 3 Application à la simulation de I'emboutissase hémisohériaue
pour des déformationsde I'ordre de 2Yo,les résultatsobtenuspar I'algorithme sont en bon

accordavecceux obtenuspar le biais du codeAbaqus (figure 3.14).
i::-
- -- ri- ----* ,- "--1--- --f
ç r
,
u.a.N.1oeezæz;
i
o
.9 0q
123456
desnoeuds
Numéro
Figure( 3.I4) : Testsur la variation de l'épaisseur pour deux élémentssoumis

à une traction axiale, comparaison avecAbaqus
décalageréapparaîtde nouveau.En conséquence,

Pour des déformationsde I'ordre de 20o/o,1e
le seul paramètrequi peut être mis en causeest la loi de comportement.Celle que nous avons
adopté ne permet pas de traiter I'aspect des grandes déformations. C'était déjà notre
hypothèsedans le chapitre2.Actuellement, un travail dans le cadre d'une thèse, est en cours
pour modéliser les lois de comportementélastoplastiqueen grandesdéformationspar des
algorithmesasymptotiquesnumériques.
3.6. Etude de la qualité de la solution
Dans cette section,nous étudionsla précisionde I'algorithmeprésenté.Pour cela nous

allons comparer l'évolution du résidu au cours du processus.Les résultatsseront comparés
avec ceux obtenus par le code Abaqus. On rappelle que, au niveau de la valeur du résidu,
Abaqus utilise une valeur par défaut de 5.10-3,si I'utilisateur ne fixe pas une valeur dansles
paramètresnumériques dans le fichier des données relatif au problème. Dans I'algorithme
présenté,la précisionest géréepar les deux paramètresde précision.Il s'agit de q et e, des
équations (2.35) et (2.36) du chapitre2. Tous les calculs présentésdans cette étude,
correspondentà des paramètresde précision suffisamment faibles 6,=sr=lQ-7 et donc aucune
opération ou phase de correction a été employée. Une étude du problème de contact en
-79-
Chapitre 3 Application à la simulation de I'emboutissaeehémisphériaue
utilisant des algorithmes asyrnptotiquesbaséesur une procédure de prédiction-conection est

en cours. On s'attend, par cette demière technique, à une augmentationde la taille des pas de
la M.A.N en utilisant des paramètresde précision moins serrésque ceux utilisés dans ce
travail.
-2.5
-3
â -1 5
J A
I
- Abaqus
- - M . A . N . 6 0 0d . d . t .
..' r { . A . N .t 5 3 6d . d . l
-----M.4.N.5766d.d.I
05tot5202530
Déplacementdu poinçon (m)
Figure (3 I 5) ' Etude de la variation du résidu,

comparaison avec le code Abaqus
Une variation du résidu au coursdu processus,peut être constatée(figure 3.15),ce qui est tout
à fait logique compte tenu du fait que l'équilibre est satisfait en tout point au cours d'un pas
de calcul de la M.A.N. Malgré cette variation, les valeurs du résidu restent en dessousde
cellesobtenuespar le code Abaqus.
-80-
Chapitre 3 Aoplication à la simulation de I'emboutissaeehémisphérioue
3.7. Conclusions
Dans ce chapitre, nous avons appliqué l'algorithme asymptotiqueprésentédans le

chapitreprécédentà la simulation tridimensionnellede I'emboutissagehémisphériqued'une
coquemince. Donc une simulation de I'emboutissageest possibleen utilisant un algorithme
basé sur la méthode asymptotiquenumérique. Le couplagede la non-linéarité géométrique,la
non-linéarité matérielle et la non-linéarité du contact, dans un cadretridimensionnel, n'a pas
altêrê l'efficacité de la M.A.N. Les résultatsprésentésmontrent clairement la précision de
I'algorithme pour prédirel'évolution des contrainteset la force de I'emboutissageau cours du
processus.L'étude de temps CPU nous a permis de localiserI'ordre optimal de troncaturedes
séries pour éviter un coût élevé du calcul du second membre. L'analyse détaillée des
paramètresde régularisationintroduits dans le modèle nous ont facilité I'ajustage de ces
paramètres pour un bon compromis entre rapidité et précision de l'algorithme. Enfin la
comparaisondes résultatsen terme de temps de calcul, avec ceux obtenuspar le code Abaqus
sont en faveur de la Méthode Asynptotique au sensde la rapidité et ceci grâce à I'aspect
automatiquede la déterminationde la taille de pas et à I'algorithme spécial du traitement du
problème de contact. Rappelons aussi que contrairement aux méthodes itératives pour
lesquellesles branchesde solutionne sont obtenuesque de manièrediscrètes(despoints), les
calculs asymptotiquesdonnent des représentationsanalytiques continues en chaque pas de
calcul.
L'algorithme présentéconserveson efficacitéquandla taille des problèmesaugmente,

ceci est très important dans la simulation des procédésde mise en forme oir les problèmes
traités sont généralementde grande taille, le nombre de points de contact est très élevé et le
raffinement du maillage est souvent nécessairepour assurer la convergenceet obtenir une
solution satisfaisante.
-81
Chapitre 4 Adoucissementde la loi de contact
Chapitre 4
Adoucissementde Ia loi de contact
Introduction 83
4.1.1.Conditions de contact avecinterpénétration 84

4.1.2.Régularisation des conditions de contact avec interpénétration 84
4.1.3.Géométriede contact 86
4.1.4.Techniquede perturbation 86
4.2. Application à la simulation de I'emboutissagehémisphérique 88
4.2.1.Gestionde la valeur d'interpénétration 88

4.2,2.Ajustage des paramètres de régularisation 89
4.2.3.Comparaisondes temps de calcul 92
4.2.4.Comparaison de la qualité de la solution 93
4.3. Conclusions 94
-82-
Chaoitre 4 Adoucissement de la loi de contact
Introduction
Dans la littérature de la simulation numérique des procédésde mise en forme des matériaux,
les conditions de contact introduisent une rigidité numérique particulière due à la
méconnaissancea priori des surfaces de contact et à la forte discontinuité des nombreuses
conditionsnon-linéairesà satisfaire.En particulieren 3D, où le nombredespoints de contact
est souvent très grand [Nagtaal et Taylor/9l]. Pour surmonter ces difficultés les codes
industriels autorisent généralementune interpénétration entre les corps rigides et les corps
déformables.
Dans la modélisation du problème de contact avec la méthode asymptotique,nous
avons vu que la loi adoptée(équation 2.15 du chapitre 2 ) n'autorisepas ce phénomène
d'interpénétration.
La dureténumériquedes conditionsde contact,dansle cadrede la M.A.N.
se manifeste par des tailles de pas de calcul relativement réduites, comparéesà celles
correspondantesaux autres non-linéarités telles que la non-linéarité géométrique et la non-
linéarité matérielle( plasticitéen déformationtotale). Ici nous précisonsque nous comparons
I'apport de la M.A.N. dans différents problèmesnon linéaires,mais comparéeaux méthodes
classiquesdu traitement du problème de contact, I'apport de la méthode asymptotique reste
très grand [Elhage-Hussein/98],[Brunelot/99],voire spectaculairesurtout lorsque le nombre
de points de contactdevienttrès grand : contact3D [Abichou et all99]. Dans ce chapitre,dans
le but d'adoucir la loi de contact, nous proposeronsune modélisation autorisant une
interpénétrationentre le corps rigide et la structure. L'objectif ici est de réduire le caractère
raide des conditionsde contact.Pour testercette nouvellemodélisation,nous I'appliqueronsà
la simulationde I'emboutissage
hémisphériqueque nous avonstraité dansle chapitre3.
4.1 Formulation du problème
La formulation du problème présentée dans le chapitre 2 reste valable dans le

développementde ce chapitre,seule la loi de contact change.Ainsi, le problème global est
décrit par l'équation d'équilibre (2.1), la loi de comportementrégularisée(2.II) et les
nouvellesconditionsde contactque nous allonsdévelopperdansce chapitre.
-83-
4.1.1.Conditions de contact avecinterpénétration
Le phénomèned'interpénétration se manifeste géométriquementpar un écart h ,entre

le corps rigide et la structure à déformer, de signe négatif. La réaction de contact dépendde
I'interpénétration(h) et le point de contactdoit satisfaireles conditionssuivantes:
o si ft20 , R=0 ( Il n'y a pas de contact) (4.1)
o si ft<O, R*0 ( Le point est en contact) (4.2)
o R.ft=O ( La condition de non-adhésion) (4.3)
Figure(4. I ) : Conditions de contact avec interpénëtration
Si l'écart ft est négatif, une réaction de contactest généréeet elle est de même sens
que la normale n, tandis que si l'écart est positif ou nul la réaction est réduite àzéro. Nous
constatonsla discontinuitéde la relation entre l'écart et la réaction du contact.Dans notre
travail, il est nécessairede détenir des relationsnon singulièrespour pouvoir appliquerune
techniquede perturbation.
4.1.2.Régularisation des conditions de contact avec interpénétration
Les conditions de contact avec interpénétration(4.1), (4.2) et (4.3) ne sont pas

analytiquesdonc elle ne sont pas développablesen sériesentières.Pour pouvoir appliquerune
-84-
technique as)anptotiquenous proposonsde régulariserces équations.Afin d'assurerun cadre

quadratique du problème pour la bonne mise en æuvre de la méthode asymptotique, nous
proposonsla régularisationsuivante:
+xn ) R = r y ( 6- h )
[rn (4.4)
lR=R.n
Où R est le module de la réaction de contact, h estl'écart entre le corps rigide et la

structuredéformable, K est un paramètrepositif (rigidité de I'obstacle)et suffisamment grand
par rapport à I'unité. â est l'écart initial entre le point de contact et la surfacerigide, n est le
vecteur normal à la surfacede contact passantpar le point de contact étudié. R est la réaction
de contact (vecteur). ry est le paramètrede régularisationde la loi de contact.
La première équationdu système(4.4) donne le module de la force de contact. Quant à

la seconde,elle déterminele senset la direction de la réactionde contactqui sont fixés par la
normale n au point de contact étudié.
Comme nous l'avions fait pour le modèle adoptédans le chapitreprécédent(section

2.1.4.b.),pour assurerune régularisationuniforme sur la surfacede contact,le paramètrede
régularisationest calculé,en chaquepoint candidatau contact,par la relation suivante:
(Rd +Khd )Rd

4,= (4.s)
6 -hd
Dans la figure (4.2), nous montrons la différence entre la loi de contact régularisée
adoptéedans les chapitresprécédents(modèle d'Elhage-Husseinet al ) et la nouvelle loi qui
autorise une interpénétration.Dans ce qui suit, pour identifier les deux modélisations nous
qualifions la loi de contact sans pénétration de 'contact dur' et le contact qui autorise une
interpénétrationde'contact adouci'.
-85-
Chaoitre 4 Adoucissement de la loi de contact
Figure (4.3) : Comparaison des lois de contact régularisées
Nous pouvons constater que pour des valeurs très grandes du paramètre K, la
nouvelle loi tend vers la loi de contact dur.
4.1.3Géométriede contact
Pour résoudre le problème de contact, les équations(4.4) doivent être complétées

(h) et n. Ces dernièressont donnéespar
grâceaux définitions des variablessupplémentaires
la géométrie du corps rigide utilisée et la cinématique du problème de contact étudié. Dans
cette partie, nous nous limiterons à un corps rigide de forme sphérique,donc nous conservons
les définitionsde ( h) et n donnéespar les équations(2.17) et (2.18)du chapitre2.
4.1.4.Techniquede perturbation
Dans cette section, nous appliquons une technique de perturbation à la loi de contact
régularisé (4.4). Nous développons en séries entières, par rapport au paramètre de
lesvariablesmisesenjeudansl'équation(4.4),àsavoir
développement<<a>>, R, ft, n, u et
2 . Après identification des développementsobtenus membre à membre par rapport à <<a >>,
nous obtenonsles systèmesd'équationssuivantes:
A I'ordre 1 :
-86-
Chaoitre 4 Adoucissementde la loi de contact
(X+uo-c-2ow , u,-2,w )
hr=
h o +P
lr-Irw -hrn,
nl =
ho* P (4.6)
n,=-fffih
R , = R o n ,+ R , no
A I'ordrep:
'
(x+no<-dw
\'-'-'o
-'!" ,t ,r-Arnl tg(or-,-2o-,w ,, -no-,-).r-,w)
'p'-p"l,l\r\-p-r'?-r"
''rL
p.r'7-t-'f
htn htn
4+p )L
4+p
ao=
(4.7)
^ KR+n^.?R+
, u^
R =-----!----!!-h -)
"p ==r:11'
2&+Kh, ft24+X4'v-,
p-l
&=4oo*4no+$no-,
r=1
Ainsi la relationnon-linéaire(a.$ est transforméeen p relationslinéaires.La linéaritéde ces
résolutionsen permetune résolutionsimple.

Nous rappelons que ces relations sont définies en chaque point de la surface
susceptibled'entrer en contact avec la surfacerigide sphérique.La régularisationdu problème
de contact permet de supposerconnus a priori les points candidatsau contact et d'imposer en
chacun de ces points un effort de contact régularisé. La régularisation des conditions de
contactrend automatiquel'étude de l'évolution du contact,c'est à dire I'analysedespositions
relativesdespoints, la prise en compte desconditionsde contactdansl'équilibre du problème
et le décollementdespoints.
Après une discrétisationpar élémentsfinis, nous obtenonsune forme matricielle de la réaction

de contact, en tout point susceptibled'entrer en contact, que nous pouvons écrire sous la
forme suivante :
t",{r"l*{,;' }
l^;Fl*,"kt,\* (4.8)
-87
chapitre 4 Adoucissementde la loi de contact
Les termesimpliqués dans cette relation sont définis dans le chapitre2. Évidemment,leurs
expressionsont changé car nous avons changé la loi de contact. Les termes de la matrice
[K,.] vont être pris en compte dans l'évaluationde la matrice de rigidité élémentaire.De
mêmepour les vecteurs{t'} ., le vecteur Fit'} 0". l'on superpose,

respectivement,
au
vecteur force élémentaireet au vecteur secondmembre élémentaire.
4.2. Application à la simulation de I'emboutissagehémisphérique
A travers la simulation de I'emboutissagehémisphérique,nous validons le nouvel

algorithme de contact. Pour cela nous reprenons le problème présentédans le chapitre 3 qui
consisteà déformerune tôle par contactavec un poinçon sphériquesoumis à un déplacement
vertical. La structure étudiée est une coque élastoplastiquede forme circulaire. Nous
conservons,ici, les mêmesdonnéesdu problème,voir la figure (3.1) du chapitre3.
4.2.1.Gestionde la valeur d'interpénétration
Dans ce modèle, nous devonsbien sûr trouver une stratégiede gestion de la valeur
d'interpénétration entre la surface rigide et la surface déformable. Dans les codes industriels
la valeur de l'interpénétration Hcrit(Allowable penetration) est généralementfixée par
l'utilisateur, sinon elle est prise par défaut. La valeur de Hcrit est calculéeen fonction de la
taille des élémentsfinis utilisés dans le problème.Dans I'algorithmeproposé,pour respecter
la valeur de Hcrit, nous proposonsde vérifier après chaque pas de calcul la condition
suivante:
R>-K.Hcrit (4.e)
Où R le module de la réaction du contact au point étudié et K le paramètreintroduit

dans la nouvelle régularisationdu contact.Après un pas de calcul si cette condition n'est pas
satisfaite,on augmentela valeur du paramètreK .
-88-
4.2.2.Aiustage desparamètresde régularisation
Les conclusions du chapitre 3 concemantle choix de I'ordre de troncature des séries et le

choix des paramètres de régularisation sont prises en compte dans ce chapitre. Ainsi les
calculsprésentésici correspondrontà un ordre 13 de troncaturede séries.
La figure (4.5) montre l'évolution de la réaction en un point de contact en fonction de
l'écart h, qui peut être négatif dans cette nouvelle régularisationdes conditions de contact.
Cette figure met en évidenceI'influence du paramètrede régularisationde la loi de contact,
comme dans le cas de la loi de contact dur. Dans cette régularisation,nous avons choisi une
valeur de 102pour le paramètreK qui gèrela valeur de I'interpénétration.
t40
t20 \
Rd,,',-0._5
100 \
Rd:1
Itr
- Rd:10
80 a\ ----Rd:5
60
\
40 \
\
20
0
0.5
h
Figure (4.4) Evolution de la réaction du contact en un point enfonction de l'écart h :

influence tle la régularisaîion de ta loi tle contact, K=102
Sur la figure (4.5), la valeurde K vaut 103, ce qui réduit remarquablement

la valeur de
I'interpénétration.
-89-
350
300
250 Itd:0.-5
- -- Rd:1
200 \ Rd--5
Rd:r0
= 150 r{d-20
100
50
0
-0.5 00.5 I 1.5 2
h
Figure (4.5) Evolution de la réaction du contact en un point enfonction de l'écart h :
influence de la régularisation de la toi de contact, K=l03
Pour comparer la nouvelle modélisation avec celle adoptéedans le chapitre précédent,nous
avons choisi une valeur de K=103 pour que I'interpénétrationne soit pas excessive(figure
4.6.).
En ce qui concerneles paramètrede régularisationRd et hd, nous avons les valeurs utilisées
dans le chapitre précédent(hd=0.05 et Rd=0.5) pour assurerune bonne estimation des
grandeurslocales.
l-----f-- I
I
...1
j
50
'
:
.:, l
ilRil 40 l..t:i clurc I
l: -.1
!
-l
- - Loi adoucie
30
t:
I
IJ
1
ll 'i
20 ll ri , .,
l;
r!
ii
10 I
':l
il"
Figure (4.6) Evolution clela réaction au ro)ro", en un point enfonction de l'écart h :

comparaison des deux lois de contact
-90-
Sur la figure (4.6), nous montrons une comparaisonentre les deux modèlespour l'évolution
de la réactionen un point de contacten fonction del'écart h .
Maintenantnous allons essayerde voir I'impact de l'adoucissementde la loi de contactsur
les résultatsobtenus.
--,a-.... Loi dure

z - - Loi douce
t
\ l
(.}. a\
€ \.
c)
(.)
L
0 0.2 0.4 0.6 0.8 I l.z

À
Figure (4.7) Evolution de laforce du poinçon enfonction du paramètre de charge :
comparaison des deux lois de contact
La figure (4.7) montreune comparaisonde l'évolution de I'effort du le poinçon en fonction

du paramètrede charge )" . Les deux algorithmes(loi de contactdur et loi de contact adoucie)
semblentdonnerles mêmesrésultats.La nouvellemodélisationn'a pas altéréla précisionde
I'algorithme dansI'estimationde I'effort global de contactou I'effort de I'emboutissage.
-91
800
8. 7rn
ç 600
5
€ soo
o
9G 4 m '--
.È
ff 3oo
o il
il
F2æ
L
LI,
il
È
ti
,
E loo i-l i ,
ti
!i
0 .r-...-..-.
...i................-1.-
0 A1
À
Figure (4.8) Evolution de la contrainte équivalentede Mises enfonction

du paramètre de charge : comparaison des deux lois de contact
La figure (4.8) illustre l'évolution de la contrainteéquivalentede Mises en fonction du

paramètre de charge ), . Contrairementà la force du poinçon qui est une force globale, la
contrainte de Mises, en tant que variable locale, est légèrement affectée par le phénomène
d'interpénétration. Un compromis doit donc être fait entre la précision et la rapidité de
I'algorithme. En effet, un adoucissementde la loi de contact permet un gain du temps de
calcul (section suivante)tout en perdantlégèrementen terme de précision.
4.2.3Comparaisondes temps de calcul
Dans cette section,nous comparonsle temps de calcul nécessairepour un déplacementdu

poinçon d'une valeur de 27.8 mm. Pour les tailles du problèmeprésentéesdans le chapitre3,
nous comparonsle nombre de pas faits par l'algorithme basésur une modélisationdu contact
dur à cellesdonnéspar le nouvel algorithme.
M.A.N
n.d.d.l
avant maintenant
600 68 39
1537 58 48
5766 55 46
'inversions
Tableau (4.I ) : Comparaison en terme de nombre d de matrices
entre les deux modèles de la loi de contact
-92 -
Le tableau (4.1) rassemblele nombre de décompositionsde matrices tangentesnécessaire

pour résoudre le problème discrétiséen 600, 1537 et 5766 degrésde liberté. Ce tableau
illustre l'apport de I'adoucissementde la loi de contact.Cet apport estplus important dansun
problèmediscrétiséen 600 d.d.l., où nouspouvonsconstaterun gain de presque50%. Lorsque
nous augmentonsla densitédu maillage,I'apport de l'adoucissementest moins important.En
fait la nouvelle loi de contact est adouci deux fois: une fois par la régularisationhlperbolique
et une fois par le coefficient de pénalisationK. Or on a vu au chapitre3, que I'adoucissement
dû à la régularisation hyperbolique est un peu plus fort avec le maillage fin, par ce qu'on a
travaillé avec un contact nodal. On peut donc comprendre que I'adoucissementdû à la
pénalisation K est plus important quand il s'applique à un contact encore "dur". C'est
pourquoi la pénalisationpermet de réduire davantagele nombre de pas avec le maillage
grossierqu'avecles maillagesplus fins. Dans la plupart des codesde simulationdu problème
de contact, le raffinement du maillage est souvent nécessairepour éviter des interpénétrations
excessiveset assurerla conversencede la solution.
4,2.4.Comparaison de la qualité de la solution
-2.5
Loi dure
-J -- Loi adoucie
â
E
'!2
-3.s
S-4 ;-.-_
F.l
-4.5
-5
-5.5
N
w .?a îwA. I
.
îvÂ
.u N
v .P
o I
À
Figure (4.9) Evolution du résidu de l'équilibre enfonction du paramètre
clecharge : comparaison des deux lois
La figure (4.9) illustre une comparaisonde l'évolution des résidusde l'équilibre donnéspar
les deux algorithmes en fonction du paramètrede charge 2. Nous constatonsque nous
pouvons réduire I'effet de la forte non-linéarité de la loi de contact tout en assurant une
solution aussipréciseque celle obtenueavecune loi de contactdure.
- 9 3-
4.3. Conclusions
Ainsi, nous avonsprésentéune nouvellerégularisationdesconditionsde contact.Cette

régularisation, au contraire de celle adoptée dans les chapitres 2 et 3, autorise une
interpénétrationentre le corps rigide et la structuredéformable.Le but est de réduire I'effet du
caractèreraide des conditions de contact. Nous avons appliqué l'algorithme obtenu à la
simulation de I'emboutissageque nous avons traité auparavanten utilisant une loi de contact
dure. Les résultatsobtenusmontrentd'abordqu'il estpossiblede traiter le problèmede contact
avec interpénétrationavec une méthode asymptotique.Ils soulignent aussi, clairement, I'effet
de l'adoucissementde la loi de contact sur le temps de calcul. Nous avons pu, grâceà cette
nouvelle régularisation,réduire le temps de calcul tout en assurantdes résultatsaussiprécis
que ceux obtenusen utilisant un algorithme de contact dur.
Nous avons constatéaussi que grâce à cette nouvelle régularisation,nous pouvons,

pour un choix spécifiquedesparamètresintroduits,modéliserà la fois la loi de contactdur" et
la loi de contactadouci".
-94-
Chapitre 5 Généralisationdes surfacesde contact
CHAPITRE 5
Généralisation de surfaces de contact
Introduction 96
5.1. Méthodes utiliséespour la représentation de surfacesde contact 96

5.1.1.Maillage des surfacesen élémentsà géométriesimple 96
5.1.2.Représentationanalytique des surfaces 96
5.1.3.Carreaux de Bézier 97
5.1.3.a Approximation polynomiale dtune courbe 97
5.3.2.b. Approximation polynomiale d' une surfoce 98
5.2. Représentation des surfacespar les fractions rationnelles 99
5.3. Représentationsdes surfacesdans le cadre de la M.A.N. 101

5.3.1.Régularisation des surfacescomposéesde surfacessimples 101
5.3.2.Traitements des surfacesarbitraires t04
5.4. Conclusions 106
95
Chapitre 5 Généralisationdes suffaces de contact
Introduction
Au cours des chapitres précédents,nous avons présenté des algorithmes qui traitent de
problèmesen présencedes conditionsde contact.Dans ces algorithmes,la géométriedes surfaces
de contactmis en jeu, étaient limités à des formes simples.Dans ce chapitre,nous nous sommes
attachés à étendre ces algorithmes au traitement des surfaces de formes géométriquement
complexes. En premier lieu, nous présentons les différentes méthodes utilisées pour la
représentationdes surfacesde contact. Dans la deuxième partie, nous montrerons comment, nous
pouvons représenterles surfacesde contact dans le cadre de la M.A.N., tout en respectantles
spécificités de la méthode à savoir la continuité des fonctions employéeset le cadre quadratiquedu
problèmeglobal.
5.1. Méthodes utilisées pour la représentation des surfaces de contact
La modélisation des outils, généralementsupposésparfaitementrigides, constitue une tâche

importanteet délicatedans la modélisationnumériquede I'emboutissage
et les procédésde mise en
forme en général. Dans la littérature, il existe trois méthodes pour modéliser numériquement les
outils, qui seront bien évidemment représentéspar leur surfacesactives.Cette section est consacrée
à la présentationde cesméthodeset à l'étudedesavantageset deslimitationsde chaqueméthode.
5.1.1.Maillage des surfaces en élémentsà géométriesimple
Pour surmonterla difficulté qui réside dans la complexité des surfacesmises en jeu, ces
dernières sont maillées par des éléments finis triangulaires ou rectangulaires.Ici, nous entendons
par éléments finis des éléments de surfaces, purement géométriques.Ces éléments sont
généralementsupposésrigides. Cette méthodeprésentedeux inconvénients:
- La modélisationdes surfacescompliquéesest délicateet diffrcile à réaliser.
- La non-régularité des surfacesmaillées engendredes problèmesnumériquesdans la résolution
du problème de contact surtout quand il s'agit de surfacesreprésentantdesrayons de congé.
5.1.2.Représentation analytique des surfacesde contact
Cette méthodeest utilisée dans la simulation des cas académiquesde mise en forme, où la
validation des logiciels de calcul a nécessitéle traitement des exemplessimples. Ceci permet
-96-
Chapitre 5 Généralisation des surfacesde contact
d'analyser les résultats expérimentaux. Sachantque cette opération est pratiquement impossible à
réaliser dans un cas industriel. Les exemples d'outillage traités sont constitués par une matrice
torique, un serre-flan en forme de couronne généralementplane et un poinçon sphérique, ou
conique à fond plat dans certainscas. Avec ce type d'outillage,où la premièreméthodeest facile
d'utilisation, on a simulé sans difficulté particulière I'emboutissagehémisphériqued'une tôle
circulaire (chapitre3 et chapitre 4). La méthode de représentationanallique des surfacesprésente
les avantagessuivants:
- Les surfacesobtenuessont de formes régulières,
- La méthodenécessiteune place mémoire dérisoire.En effet, seul le stockagede quelques
paramètres(rayons,centre,...) est nécessaire,contrairementà la premièreméthodequi nécessitele

stockagedescoordonnéesde tous les sommetsdeséléments,ainsi que leursconnectivités.
- L'outillageesttrès facile à gérer.
La limitation de cette méthode réside dans le fait que les surfacesrencontréesdans les problèmes
industriels n'admettentpas touj ours des représentationsanalyiques.
5.1.3.Carreaux de Bézier
Faute d'avoir une représentationanalytique ( une carte ) exacte ou approchéedes outils, les
méthodesd'approximationsde surfacesse sont imposées.lnspiré du métier de I'emboutisseurqui
utilise une méthode très particulière pour gérer les courbes ( qui ne sont que des sections de
surfaces),P. Bezier a mis en æuvre dans un premier temps, une méthode d'approximationbien
adaptéeà ce métier. Cette méthode a été étendueensuitepour traiter le cas des surfaces[Farrin/90].
5.1.3.a. Approximation polyn omiale d' une courbe
Etant donnée une courbe C, nous pouvons approximer les coordonnées d'un point
appartenantà cette courbepar:
N
x(r)=f d,t'
; -|
N
yu)=L Ft', (s.l)
i=l
lr'
zQ):Z r J'
-97 -
Où I'on a introduit un paramètre re[O,t], où ti représentet àla puissance i, N le degré du
polynôme,a, , f , ,2/, sescoefficients

vérifiantx(O)=ao,y(0)=Fo,z(0)=yo les coefficients
de
I'origine de la courbe. Les coordonnées de I'extrémité de la courbe vérifient
x(l):st , y(l)= Ft z(O)=yr,.

Si on note P(r) le vecteurdescoordonnéesd'un point quelconquede la courbe,nous avons:
('(0')
P(1)1y@ |
\'(t))
Cequi peutsemettresousla forme:
P ( t ) : b i t i , i = | . . . N ,b , = b , , p , , y , f
Cetteéquationdéfinitparfaitement par (l).
la courbeparamétrée
5.1.3.b. Approximation polynomiale d' une surface
En partant du principe q'une surface peut être obtenue en faisant varier les points d'une
courbe dans un domaine bien défini. Chaque point de cette courbe va décrire une courbe
appartenant à la surface. Si la première, dite courbe génératrice,peut être approximée par un
polynôme de degré N, chaque courbe engendréepar le mouvement de la génératrice,peut être
approximéepar un polynômeM. Les deux courbesdécritespar le mouvementdes points extrêmes
de la génératricesont appeléescourbesdirectrices.Ces courbessont paramétréespar une variable
mefÙ,ll. Bien évidemmentpour m=0 on retrouvela génératriceinitiale et pour m=I on retrouvela
génératricefinale. La surface ,S est paramétréepar t et m.Un point appartenantà la surfaceest

notê P(t,m), il s'exprimepar:
NM
rQ,m1=llB,,r'^t (s.2)
i=t j=r
Une telle représentation répond bien à notre cas car elle s'adapte bien au développement
asymptotique. Elle permet à partir d'un ensemble de points donnés de définir une surface
quelconquequi approcheles surfacesréelles de I'outil. Ces coefficients ô'' et Bu sont calculésà
partir des coordonnéesde points choisis sur la surface réelle. Ceci permet d'avoir toutes les
informations utiles en chaquepoint de la surface (le plan tangent et la normale à la surfaceen un
point).
-98-
Chapitre 5 Généralisation des sur.facesde contact
5.2. Représentationdes surfacespar les fractions rationnelles
Si nous résumonsce que nous avons présentédans la sectionprécédente,nous pouvons dire

que dans la modélisation géométrique pour la représentationdes surfaces et des courbes, deux
conceptsmathématiquesnon linéairessont généralementutilisés:
. La modélisation implicite qui permet de représenterfacilement les quadriques simples (
ligne, cercle,conique,plan, sphère,...)mais ne permet absolumentpas la représentationde formes

quelconquesqui peuventêtreutiliséesdansla conceptionou la fabrication.
. Les polynômes paramétriques (splines)représententbien les courbesquelconquesmais
sont relativement imprécis pour la représentationdes quadriquessimples.
Pour éviter les inconvénientsdes deux approchesci dessus,nous proposonsd'utiliser les B-

splinesrationnelles et non uniformes sousforme canoniquepour représenterla géométrie.
Pour générerdesprofils, trois typesde courbepeuventêtreutilisés:

- Un segmentde droite.
- Un arc de cercle.
- Une splinequelconque.
cesgéométriespeuventêtredéfiniespar la formule suivante[Tiller/83]:
Mathématiquement,
ln,,o1t1.t',.q
C(t)=(x(t),y(l),O;=-=1-- ( 5.5)
ln''oçt1'n'
où .Bsont les points bidimensionnelsappeléspoints de contrôle,t est le paramètre( a9<b ) aveca
et b fixes tel que : 03a<b
-B;,rsont les coefficientsdes polynômesde la variable t d'ordre fr et de degrét-1. Ils sont appelés
les fonctionsde base.
Ces fonctionssont complètementdéfiniespar I'ordre k etle vecteur {rr};]f ori '
a- f t z-...- t,,4,,*r1t,,* 11...

tr9t n*r=...=tn*k:b .
S'il existeun nombrepositif d tel que
-99-
Chapitre 5 Généralisationdes surfaces de contact
t,*r_t,=d polut k<<n

rl
le vecteur f, Iest alors uniforme.
La valeur de Bi,rpour / fixé est calculéen utilisant I'algorithmerécursifde Cox-De Boor.
homogènes( h,> )pour tout i.

ftisont les coordonnées
Les surfacessont représentéesmathématiquementpar :
Ilq,-t'l B,,t?)hrPu
S(l)=(x(r,s),y(t,s),2(t,s))=# (5.6)
Ii=r la,-(")Bj,r(t)htj
j=r
Avec R1 un ensemble de nxm points de contrôle tridimensionnel, hii un ensemble de nxm

coordonnéeshomogènes.
Dans le cas des splines quelconques,en plus de la position des points spécifiés,
I'utilisateur a quelqueslatitudesquant au nombre de points et à I'ordre des B-splines devant être
utilisés. L'ordre détermineles propriétésde continuitéde la courbeet par la même sa forme. Si les
splines sont utilisées, leurs dérivés jusqu'à k-2 sont continues.Par exemple si fr=4 (B-spline
cubique), la tangente et la courbure sont continues, par contre si k=3 (quadratique) seule la
tangenteest continue.De façon généraleplus I'ordre est élevé,plus lourd est le systèmeà résoudre.
Aussi I'utilisateur doit toujours utiliser l'ordre minimum pour obtenir le niveau correct pour la
continuité.
Le nombre et la position des points spécifiés par I'utilisateur sont fondamentaux pour la
détermination de la forme de la courbe. Plus de points donnent un meilleur contrôle de la forme
mais conduisentà desvolumesde calculsplus importantes.
Dans le cadre de la Méthode Asymptotique Numérique, la représentationdes surfaces
de contact au moyen de ce processuss'avèretrès intéressantedans le sensor) elle permet de traiter à
la fois les surfacesimplicites et les surfacesarbitraires. Par ailleurs, ce genre de représentations
s'adaptebien au développementsasymptotiques.Aussi nous sommes convaincus,d'après nos
expériencesavec les séries polynomiales et fractionnelles, que ces derniers sont plus efficaces et
plus rapides.
-100-
Chapitre 5 Généralisation des surfaces de contact
5.3. Représentationdes surfaces de contact dans le cadre de la M.A.N.
Introduction
La modélisation du problème de contact dans le cadre de la méthode asymptotique

numériquefait intervenir la géométriedes surfacesrigides (équations2.17 et 2.18 du chapitre2).
Cette caractéristiquen'a pas générêde difficultés particulières dans le cas de géométriessimples
tellesque les lignes,cercles,plans et sphèresfElhage-Hussein
e/ al./98]1,fBrunelot/99]et [Abichou
et al. /01]. Ceci n'est pas le cas lorsqu'il s'agit des corps de géométriescomplexes où les
représentationsanalytiquesne sont pas toujours évidentes.Il suffit de prendre comme exemple la
simulation des procédésde mise en forme dansle domainede I'automobile.Une représentationde
ces outils nécessitedonc la mise en jeu de différentesfonctionscapablesde décrireces différentes
formes.
5.3.1.Régularisation des surfaces de contact.
Le problèmequi se pose au niveau de la simulationdu problèmede contact,en présencedes

surfaces de formes complexes, est le passaged'un point de contact d'une surface à une autre.
Autrement dit au niveau de la connexion ou de la jonction de surfaces différentes. Comme le
montre la figure 5.1, un point de contactoccupantla position i se trouve aprèsun certaintempsà la
position i'en quittant la surface de contact (figure 5.1.a) ou en pénétrantla surface de contact
(figure5.1.b)
i
(a)
Figure (5.1) Correction de la position d'un point de contact par rapport à la sudace rigide:
(a):cas de décollement, (b) cas d'interpénëtration
Kt--"tT-rà
-101- ç/
Dans les méthodes itératives de traitementsdes problèmes de contact, la phase de correction,

permet de rectifier la position de ce point en le projetant sur la surfacede contact (position i" figure
5.1). De toute évidence,cettetechniquegénèreun tempsde calcul prohibitif, si I'on se met dansle
cas d'un problème de mise en forme où le nombre des points de contact est très grand et le
changementd'allure des surfacesde contactest très fréquent.
Dans le cadre de la MAN, cette procédure de correction ne peut être adoptée.Et nous ne
pouvons pas arrêter le calcul chaque fois qu'un point de contact quitte une surface pour entrer en
contact avec une autre, ceci va pénaliserla méthode.Des travaux sont en cours pour I'application de
la MAN avecune techniquede prédictioncorrectionaux problèmesde contact.
Pour surmonterces difficultés, sanspasserpar une procédurede correction,nous proposons
de procéder par une technique de régularisation des fonctions qui décrivent les surfaces,dont la
géométrieprésentedes discontinuités.
L'idée de base dans cette technique est de construire les surfacescomplexes à partir de
surfacessimples.Ainsi une surfacede forme compliquée va être construiteà partir d'une succession
de raccordement de surfaces simples. Nous prenons donc la surface compliquée, nous la
subdivisons en surfacessimples; le passaged'une surfaceà I'autre est assurépar la régularisation.Si
nous nous plaçons à une échelle bidimensionnelle,toute surfacecompliquée sera construite à partir
de raccordementsligne - ligne ou ligne - arc de cercle. Il est évident que cette tache devient plus
difficile si nous nous plaçonsdans un cadretridimensionnel.
L'autre motivation d'utiliser cette procédurede régularisationest la limitation des courbesde
Bézier pour représenterdes surfacesqui présententdes rayons de congé.Dans les procédésde mise
en forme, ces rayons de congé ont des dimensions préciseset constituentdes zones où la structure
qui va entrer en contact a un comportementtrès sensiblesurtout lorsqu'il s'agit de congésde rayons
très faibles.
Aux chapitres3 et 4, nous avonsprésentédes algorithmesde simulation de I'emboutissage
hémisphérique.Le traitement du problème de contact était limité à une surfacerigide qui représente
le poinçon. Ici nous nous sommesattachésà introduire une surfacerigide qui représentela matrice
dans le procédé de I'emboutissage.Cette matrice est constituée d'une partie plane et d'une partie
arrondiedéfinie par un rayonde congé r (figure 5.2).
-102-
Géométrie: Propriétés
Mécaniques
:
w
R:59.18mm E:206.8Gpa
P=50.8mm v=0.3
e:1mm ov=170MPa
n:4.75
Figure(S.2): simulation en présence d'une matrice evec un rayon de congé
Pour représentercette matrice, nous proposons de la construire à partir d'un ensemblede

raccordsligne-cercle.Nous définissons,dansun repèrelié au raccord,tout point M du raccordpar
(x)
tl
le vecteurMl Y I. Géométriquement
nousproposons
de représenter par la
un raccordligne-cercle
tl
lz)
fonction suivante:
-x2 ,2
Z_ + Î',I-
'X-r ( 5.3)
où r est le rayon de congé et q un paramètrede régularisationde la géométrie,il est positif et
suffisammentpetit.
Dansuneconfiguration
quelconque,
un pointestdéfinipar le vecteurrn[;l tel que:
l')
-=Fii).(,) (s.4)
Où [R]représente la matrice de rotation du repère lié au raccord par rapport au repère global. Le
vecteur (t) représentela translation du repèrelié au raccord par rapport au repèreglobal.
-103-
Chapitre 5 Gënéralisation des suffaces de contact
La régularisationproposéea été testéepour s'assurerqu'elle permet bien une représentation

d'un raccord ligne-cercle. Par contre son implémentation dans nos algorithmes de contact sera
réaliséeultérieurement.
5.3.2Traitementsdessurfacesarbitraires
Pour faciliter la tâche, nous allons nous placer dans un cadre 2D. Nous considéronsune
surface de contact représentéepar une courbe de Bézier d'ordre n. Certes I'utilisation de la
représentationcomplète de la surface va générer un temps de calcul énorme. Dans le cadre de la
MAN, pour minimiser les développementsdes calculset le temps de calcul, nous proposonsI'idée
suivante.
En se basant sur le fait que pour chaque point de contact, on construit une matrice de rigidité de
contact [f *] , nous n'allons plus définir la surfacede contactune fois pour toute au début du calcul
mais nous allons définir pour chaque point de contact une surface élémentairede contact. Cette
surfaceélémentaireest définie ,elle aussi,par une courbede Bézier mais d'ordrem limité (*<<p).
Cette surfaceélémentaireest un morceaude la surfaceglobale.C'est le morceaule plus proche du
point de contact étudié. Après un pas de calcul, la surface élémentairede contact est mise à jour.
Ceci ne changerien en terme de temps de calcul puisque la matrice [Kr.] est évaluéeau début de
chaquepas.
Exemple:
,:
rigide
Figure 5.3 représentation de surfaces de contact définiespar dessurfaces

arbitraires d'ordre élevé
- 104-
Chapitre 5 Généralisation des surfaces de contacl
Ainsi, si I'on choisit de limiter la surfaceélémentaireà I'ordre 7 par exemple,le point de

contact {, sera étudié par rapport à la surface élémentaireconstruite par le polygone défini par les
points 4, Pr, P3, P4, P, Pu et Pr. Quant au point lr,ll seracalculépar'rapportà la surface
construitepar le polygonedéfini par les points de contrôle P, Po,P, , Pu, P, , Pret Pn.Le point le
plus procheest 'au milieu'de la liste des points,qui définissentla surfaceélémentaire.Ceci permet
d'assurerque cette surface élémentaire représenteefficacement la surface de contact, lorsque le
point étudiése déplaceen amontou en aval du point le plus proche.
Nous pouvons constaterque si on suit le trajet d'un point de contact,nous trouveronsque ce trajet
est guidé par un ensemble de surfaces élémentairesqui sont en recouvrement entre elles. Ces
surfacesen recouvrementreprésented'une façon assezprécise la surface globale, si on choisit un
ordre suffisant pour la surfaceélémentaire.
Alsorithme:
Au débutde chaquepas de calcul et pour chaquepoint de contact:

- On cherche dans la liste des points du polygone qui définie la surface globale, le
point le plus proche.

- On construit la surfaceélémentairecorrespondante.
- On construitla matrice [K,"] calculéepar rapportà cettesurface
Cette technique que nous venons de présenterpeut être appliquéede la même manière pour
par les fractionsrationnelles(section5.2). Ceci va nous permettrede traiter à la
les représentations
fois les surfacesde formesquelconqueset les surfacesde formessimplesou implicites.
-105-
5.4.Conclusions
Les procédés de mise en forme des matériaux font intervenir des outils de géométrie
généralementcomplexe. Au cours des simulations numériquesde ces procédés,ces outils sont
souvent représentéspar leurs surfacesactives.La représentationde ces surfacesconstitue une tâche
très délicate. Une représentationdoit être assezprécise pour décrire véritablement la forme des
outils et par la suite la conceptiondu produit qu'on désireobtenir à la fin de I'opérationde la mise
en forme. La représentantgéométriqueest un facteur déterminant de I'efficacité de I'algorithme de
contact.
Dans le cadrede la MAN, la nécessitéd'avoir des représentations
analytique,donc développables
en séries entièresnous a poussé à proposer des techniquesspécifiquespour le traitement des
surfacesde formes complexes.Nous avons proposé,d'abord une méthodede régularisationqui
consiste à construire les surfaces compliquées par un ensemble de raccordementsde surfaces
simples. Dans le cas des problèmes qui utilisent des surfaces arbitraires: surfaces de Bezier
(Polynômes) ou surfacesB-Surfaces ( fractions rationnelles), nous avons à traiter ces surfaces,qui
ont des ordres élevés,en les décomposanten surfacesd'ordres limités ou moins élevéspour réduire
le nombre de variables intermédiaires et par la suite assurer une bonne mise en ceuvre de la
Méthode Asymptotique.Les théoriesprésentéesdans ce chapitrene sontpas encoreimplémentées
dans les algorithmesde contactprésentésdans les chapitresprécédents.Ceci constitueraI'objet de
nos prochainstravaux.
-106-
Conclusions qénérales et perspectives
Conclusions généraleset perspectives
Dans ce travail de thèse, nous avons présenté des algorithmes asymptotiques

numériquespour la simulation de l'emboutissageà froid. Une simulation de I'emboutissage
en 3D est, donc possible, en utilisant un solveur asymptotiquenumérique.On pouvait se
demander si I'efficacité de la M.A.N. était altérée par le couplage des non linéarités
géométriques,matérielles et de celles dues au contact unilatéral : d'après les résultats obtenus
dans cette thèse,cette efficacité est préservée.Pour réaliserdes calcul de sériesentières,il a
été nécessairede régulariser la loi de contact unilatéral: il est clairement établi que cette
régularisation n'affecte pas du tout la prédiction des quantités globales ( force
d'emboutissage)et assezpeu celle des quantitéslocalescomme les contraintes.La méthode
asymptotiques'est montrée bien plus performanteque le code commercialAbaqus, qui est
réputé être un bon code de calcul implicite d'emboutissage.La comparaisondes résultatsen
terme de temps de calcul, avec ceux obtenuspar le code Abaqus sont en faveur de la Méthode
Asymptotique au sensde la rapidité. Ceci est obtenu grâce à la gestion automatiquede la taille
de pas et à un traitement spécifique du contactunilatérale. Rappelonsaussi que contrairement
au méthodesitératives pour lesquellesles branchesde solution ne sont obtenuesque de
manière discrète (des points), les calculs asymptotiques donnent des représentations
analyiques continuesà chaquepas de calcul.
L'algorithme présentéconserveson efficacité quand la taille des problèmesaugmente:

ceci est très important dans la simulation des procédésde mise en forme où les problèmes
traités sont généralementde grande taille, le nombre de points de contact est très élevé et le
raffinement du maillage est souvent nécessairepour assurer la convergenceet obtenir une
solution satisfaisante.
Nous avons proposé aussi une nouvelle régularisation de la loi de contact que nous
avons qualifié de régularisation adoucie.Avec cette nouvelle modélisation nous avons montré
que nous pouvons traiter des problèmes de contact avec une possibilité d'interpénétrations
entre les corps en contact. Le nouvel algorithme a été testé sur un problème de simulation de
I'emboutissagehémisphérique.Les résultatsobtenusmontrent bien que nous pouvons réduire
I'effet du caractèreraide des conditions de contact, tout en assurantune bonne précision de
I'algorithme et surtout en réduisant le temps de calcul. Nous avons constaté aussi que cette
nouvelle régularisationnous permet,pour un choix spécifiquedesparamètres,de modéliserà
la fois la loi de "contactdur" et la loi de "contactadouci".
r07-
Conclusionsgénérales et perspectives
Dans le but de simuler le procédéde I'emboutissage

et les procédésde mise en forme
en général, en tenant compte des différents outils qui peuvent être mis en jeu, à savoir les
serre-flans, les matrices, les joncs de retenue, nous avons présentéles différentes méthodes
que nous pouvons utiliser dans le cadre de la MAN pour représenterles surfacesde formes
compliquées.Ces méthodesprennenten compte les spécificitésdesMéthodesAsymptotiques,
à savoir la continuité des fonctions et des relations employées et 1e cadre quadratique du
problème étudié pour assurerune bonne mise en æuvredes techniquesde perturbation.
Dans notre équipe, on vise prochainement le perfectionnement des algorithmes

présentésen développant des lois de contact entre des corps déformables et en introduisant
des lois de comportement plus réalistes, qui prennent en compte l'aspect des grandes
déformations et le phénomènede la déchargeélastique.Pour élaborerun logiciel prototlpe de
simulation de I'emboutissage
et desprocédésde mise en forme en général,nous allons inclure
dans les algorithmes obtenus I'aspect 'micro-structure' qui a un rôle important dans les
procédésde mise en forme. Nous allons aussi offrir la possibilité de traiter ces problèmesà
I'aide desprocéduresde 'sous-structuration'et descalculsparallèles.
1 0 8-
R éfé r e n c es h ih li. o qr an h i.ou es
Référencesbiblioeraphiq ues
[ABAQUS/95] ABAQUS. Theory and users manuals (version 5.3). Hibbitt, Karlsson and
Sorenson,Inc., 1080Main Street,Pawtucket,zu 02860,USA, 1995.
[ABAQUS/98-11 ABAQUS. Theory and users manuals (version5.8). Hibbitt, Karlsson and
Sorenson,Inc., 1080Main Street,Pawtucket,RI02860, USA, 1998.
[ABAQUSI9S-21ABAQUS. Examplesproblem manuals (versionJ.8). Hibbitt, Karlsson and

Sorenson,Inc., 1080Main Street,Pawtucket,RI02860, USA, 1998.
[Abichou et all99] H. Abichou, J. Brunelot, H. Zahrowi et M. Potier-Ferry,First steps

towards a simulation of material forming processes by an Asymptotic Numerical Method,
ECCM'99, Munich, Allemagne, 1999.
[Abichou et all0ll H. Abichou, H. Zaltrouni et M. Potier-Ferry, First attempt to simulate

deep-drawing processes within an Asymptotic Numerical Method, Computer Methods in
Applied Mechanicsand Engineering,2001, Soumis.
[Afzali et Bouhelier&4l M. Afzali et C. Bouhelier, Emboutissage,intérêt de la simulation et

produits existants,Recueil de conférences: Calcul en emboutissage,applications industrielles
pour la conceptiondespièces,PublicationCETIM, 1994.
[Ammar/96] S. Ammar, Méthode asymptotique numérique perturbée appliquée à la

résolution des problèmes non linéaires en grande rotation et grand déplacemenf,Thèse de
doctoratde I'universitéde Laval, Québec,1996.
[AndersonlS2lB. S. Anderson,A numerical study of the deep drawing process, Numerical

Methods of Industrial Forming Processes,(Eds) Pittman et al., Pineridgepress,Swansea,U.
K., pp.:109,1982.
[Avitzur/831Handbook of metalformingprocesses,Wiley-Intersciences,
1983.
l{zrar et all93lL. Azrar,B. Cochelin,N. Damil et M. Potier-Ferry,An asymptotic-numerical

method to compute the post-buckling behavior of elastic plates and shells, Intemational
Journalfor NumericalMethodsin Engineering,36, pp.: l25l-I277,1993.
[Baker et all96] G. A. Baker et P. Graves Morris, Padé Approximants, Encyclopedia of

Mathematicsand its applications,2ndEdition, CambridgeUniversityPress,59,1996.
[Bathe et Chaudhary/851 K. J. Bathe et A. Chaudhary,A solution methodfor planar and

axisymmetric contact problems,International Journal for Numerical Methods in Engineering,
21, pp.:65-88,1985.
[Bathe et Chaudhary/86] K. J. Bathe et A. Chaudhary,A solution methodfor static and

dynamic analysis for three dimensional contact problem with friction. Computers and
Structures,
6, 24,pp.:855-873,
1986.
-109-
Références biblio sraphiques
[Bathe et Mijailovich/88] K. J. Bathe et S. Mijailovich, Finite elementanalysis of frictional

contactproblems, Journal de MécaniqueThéoriqueet Appliquée, Specialissue, 1,7, pp.:31-
45,1988.
[Bathe/96] K. J. Bathe, Finite Element Procedures,Prentice-Hall,Engewood Cliffs, New

Jersey07632,1996.
lBatoz et all&9l J. L. Batoz, P. Duroux, Y. Q. Guo et J. M. Detraux, An fficient algorithm to

estimatethe large strains in deep-drawing,NUMIFORM'89, Balkema,pp.:383-388,1989.
lBatoz et Dhattlg0l J. L. Batoz et G. Dhatt, Modélisation des structurespar élémentsfinis,

(Eds) Hermes,Paris, 1990.
lBatoz et all9ll J.L.Batoz, Y. Q. Guo et F. Mercier, The inverseapproachincluding bending

effectsfor the analysis and design of sheet metal forming parts, NUMIFORM'95, Balkema,
pp.:661-661,1989.
[Betlet et all90] M. Bellet, E. Massoni et J. L. Chenot, Calcul des membranes

élastoplastiques et élasto-viscoplastiques,Application à I'emboutissage et au formage
superplastique, Physique et Mécanique de la Mise en Forme des Métaux, Moussey et
Franciosi,(Eds) Paris,CNRS, 1990.
[Ben Dhia/89] H. Ben Dhia, Modelisation and resolution by a penalty duality method of
uniI ateral contact prob Iem.NUMIFORM89, I 989.
[Ben Dhia/901H. Ben Dhia, Modelling and numerical approachof contactand dry friction in
simulation of sheet metal forming, WCCMZ, 2, pp.:779-782, 1990.
[Boivin/90] M. Boivin, Prévision du retour élastique lors du pliage des tôles, Physique et
Mécaniquede la Mise en Forme desMétaux, Moussy et Franciosi,(Eds)Paris,CNRS, 1990.
[Boutyour et all93l E. H. Boutyour, B. Cochelin, et M. Potier-Ferry,Calculs despoints de

bifurcation par une méthodeasymptotiquenumérique, Compte rendu du premier Congrès
National de Mécaniqueau Maroc, pp.:371-378,1993.
[Boutyour/941 E. H. Boutyour, Méthode Asymptotique Numérique pour le calcul de

bifurcations: Application aux structures élastiques,Thèsede doctorat de l'université de Metz,
t994.
[Braikat/95] B. Braikat, Méthode asymptotiquenumériqueet fortes non linéarlés, Thèsede

doctoratde la faculté des sciences,Casablanca,1995.
[Brunelot/991 J. Brunelot, Simulation de la mis en forme à chaud par la Méthode

AsymptotiqueNumériqrze,Thèsede doctoratde I'universitéde Metz,1999.
[BruneU88] M. Brunet, A finite elementmethodfor unilateral contactand friction problem

involvingfinite strain and large displacement,Joumalde MécaniqueThéoriqueet Appliquée,
Specialissue,l, 7, pp.:209-220,1988.
-110-
Ré.férences bi bI i ogr aphiques
[Brunet/901 M. Brunet, Modélisation numérique par élémentsfinis d'un emboutissage

profond, Ecole d'été d'Oléron, dirigé par F. Moussy et P. Franciosi,(Eds ) Paris, CNRS,
pp.:516-521,1990.
[BudianskylT4l B. Budiansky, Theory of buckling and post-buckling behaviour of elastic

structures,Advancesin Applied Mechanics,14, pp.'.65, 1974.
[Bussy et Vauchez/89] P. Bussy et P. Vauchez, Simulation numérique du processus

d'emboutissage,g'-" Congrès Françaisde Mécanique,Metz 5-8 septembre,pp.: 528-531,
1989.
[Biichter et Ramm/92] N. Btichter et E. Ramm, Shell theory versus degeneration- a

comparison in large rotation finite element analysis, International Journal for Numerical
Methodsin Engineenng,34, pp: 39-59, 1992.
[Cadou/97] J. M. Cadou, Méthode asymptotiquenumériquepour le calcul des branches

solutions et des instabilités dans les fluides et pour les problèmes d'interaction
fluide-structure,Thèsede doctoratde I'universitéde Metz,1997.
[Causse-Massoni187'lE. Causse-Massoni, Modélisation numériquede l'emboutissagepar

fims, Thèsede doctoratde I'Ecole NationaleSupérieuredesMines de Paris, 1987.
elements
[Chabrand et al l88l P. Chabrand,Y. Pinto et M. Raous, Numerical modelling of friction for

metal forming processes, EUROMECH, Modelling of metal forming processes, Sophia
Antipolis, pp.: 93-99, I 988.
[Charlier et Cescotto/8Sl R. Charlier et S. Cescotto,Modélisationdu phénomènede contact

unilatéral avecfrottement dans un contexte de grandes déformations,Journal de Mécanique
Théoriqueet Appliquée,Specialissue,I, 7, pp:.177-192,1988.
[Chenais et Rousselet /881 D. Chenaiset B. Rousselet,Conceptionoptimale de structure

assistéepar ordinateur,Ecole d'été,Sophia-Antipolis,1988.
[Chenot/93f J. L. Chenot, La modélisationnumérique des procédés de mise en forme, La

Revuede Métallurgie - ClT/Scienceet Génie desMatériaux,pp.:l567-1576,1993.
[Chenot/97| J. L. Chenot,New tends in finite elementmodeling of metalforming processes,

ComputationalPlasticity Fundamentaland Applications, Barcelone,CIMNE, pp.:209-223,
1997.
[Chertier/97] O. Chertier, Contact et frottement entres solides déformables en grandes

déformations,Thèsede doctoratde I'universitéd'Aix Marseille2,1997.
TCNRS/901Physiqueet Mécanique de la mise en forme des métaux, Ecole d'été d'Oléron,

dirigé par F. Moussyet P. Franciosi,(Eds ) Paris,CNRS, 1990.
[Cochelin et all94] B. Cochelin,N. Damil, et M. Potier-Ferry,Asymptoticnumericalmethods

and padé approximantsfor nonlinear elastic structures,InternationalJournal of Numerical
Methodsin Engineeing, 37, pp.: I I 87-1213,1994.
- 111-
Ré.férences bib li o sr aphiques
[Cochelin/g4|B. Cochelin.A path-following techniquevia an asymptotic-numericalmethod,

Computerand Structures,53, 5 , pp.:1181-1192,1994.
[Cochelin et Compain/99] B. Cochelin et C. Compain, Une Méthode Asymptotique

Numériquepour la dynamiquenon-linéairetransitoir€,4" " Congrèsde Calcul desStructures,
P. Ladevèzeet M. Raous,18-21mai, Giens,pp.:239-244,1999.
D. Guédra-Degeorges,
[Combescureet all89l A. Combescure,A. Millard et P. Verpeavx,Numerical methodsin the

CASTEM systemfor the treatment of the contactproblems, Journal de Mécanique Théorique
et Appliquée, Specialissue,l, 7, pp;249-261, 1988.
[Courtois et Ricard/88] P. Courtois,A. Ricard,Finite elementanalysisof industrialfriction

contactproblems, Joumal de MécaniqueThéoriqueet Appliquée,Specialissue,l,7,pp.:263'
276,1988.
[Crisfield/S1] M. A. Crisfield, A fast incrementaliteration solutionprocedure that handles

snap+hrough,Computersand Structures,13, pp.:55-62,1981.
[Crisfield/83] M. A. Crisfield, An arc-length methodincluding line search and acceleration,

lnternationalJoumalof Numerical Methods in Engineering,19, pp.:1269-1289,1983.
[Crisfield/9l] M. A. Crisfield, Non-linear finite elementanalysis of solids and structures,

John Wiley & Sons,1991.
[Damil et Potier-Ferry/90] N. Damil et M. Potier-Ferry,A new methodto computeperturbed

bifurcation: Application to the buckling of imperfect elastic structures,International Joumal
of EngineeringSciences,26, 9, pp.:943-957, 1990.
[Dhatt et TouzotlS4l G. Dhatt et G. Touzot, Une présentation de la méthode des éléments

finis, Maloine, Paris, 1984.
[Durville/93] D. Durvllle, Modélisation Mécaniqueet simulation numériquede l'avalement

dans le procédé d'emboutissage,Thèsede Doctoratde I'Ecole Centralde Paris, 1993.
[Curnier et Alart/88] A.Curnier, P. Alart, A generalized Newton method for contact

problems with friction, Journal de Mécanique Théorique et Appliquée, Special issue, l, 7,
p p . : 6 7 - 8 21, 9 8 8 .
[Duroux et Batozlg0l P. Duroux et J. L. Batoz, Une méthoded'évaluation des déformations

dans les tôles minces embouties,Physique et Mécanique de la Mise en Forme des Métaux,
Moussy et Franciosi,(Eds ) Paris , CNRS, 1990.
[Elhage-Hussein/981A. Elhage-Hussein,Modélisation des problèmes de contact par une

méthodeasymptotiquenumériqne,Thèsede doctoratde I'universitéde Metz, 1998.
[Elhage-Hussein et all98l A. Elhage-Hussein,N. Damil et M. Potier-Ferry,An Asymptotic

Numerical Algorithm for frictionless contactproblem, Revue Européennedes ElémentsFinis,
7 , p p . :I l 9 - 1 3 0 , 1 9 9 8 .
-ltz-
Références b iblio gr aPhiques
[Elhage-Hussein et all00] A. Elhage-Hussein,N. Damil et M. Potier-Ferry,A Numerical

continaation method based on Padé Approximanls, Intemational Journal of Solids and
Structures, -7001,2000.
37, pp.:6981
[El Mouatassim et Jameux/94] M. El Mouatassimet J. P. Jameux,La simulation numérique

de I'emboutissage chez RENA(ILT, validation de logiciels explicites. Recueil de
applicationsindustriellespour la conceptiondespièces.
conférences: Calcul en emboutissage,
PublicationCETIM, 1994.
[Gallagher/751 R. H. Gallagher, Perturbation procedures in nonlinear finite element

structural analysis, Computational Mechanics-Lecture Notes in Mathematics,
Springer-VerIag,Berlin, 461, pp.:15 -89, 1975.
[Farrin/90] G. Fanin, Curves and Surfacesfor ComputerAided GeometricDesign. Computer

Scienceand Scientific Computing,secondedition, 1990.
[Fourment et at/97] L. Fourment,K. Mocellin et J. L. Chenot,An implicit contoctalgorithm

CIMNE,
for 3D simulation of the forging process, Computational Plasticity, Barcelone,
-8
Fundamentalsand Application,pp.:873 I 7,1997'
du
[Franchomme et al t89l P. Franchomme,A. Ricard, J. Oudin et Y. Ravalatd,Simulation
contact 2D dans SfSfUS par la méthode du Lagrangien Augmenté, SttuCoMe 89, (Eds)
Paris,pp.: 187-198,1989.
la
[Gelin/9Qf J. C. Gelin, Analysede problèmesen mise enforme, Physiqueet Mécaniquede
Mise en Forme desMétaux, Mousseyet Franciosi,(Eds) Paris,CNRS,pp.:359' 1990.
[Gelin et all9|l J . C. Gelin, L. Boulman, et P. Boisse,Quasi-staticimplicit transientanalyses

of sheet metal forming using a C0 three node shell element, Journal of Material Processing
Technology,50, l-40, pp.: 54-69,1995.
[Germain/8s| Y. Germain,Modélisationpar éléments finis d'écoulementviscoplastiqueavec

Nationale des
frottement - Application au forgeage à chaud. Thèse de doctorat de I'Ecole
Mines de Paris,l985.
[Gosh et Hecker/7s| A. Gosh et S. S. Hecker, Failure in thin sheetsstretchedover rigid

punches,Metallurgical Transactions,6A, pp.: 1065-1074, 1975.
[Gotoh et Ishise/78| M. Gotoh et F. Ishise, A finite element analysis of rigid-plastic

deformation of theflange in a deep drawing process based on fourth-degree yield function I,
lnternationalJoumalof MechanicalSciences,20, pp.:423,I978.
[GotohiSQlM. Gotoh,A finite elementanalysisof rigid-plastic deformationof theflange in a

deep drawing process based on fourth-degree yield function II, International Journal of
MechanicalSciences,22, pp.: 367,1980.
[Grolade et Cristescu/83| D. Grolade,C. Cristescu,Résolutiondu problème de contact entre

solidesélastiques.CETIM, Senlis,1983.
113
Ré-férences bib lio gr aphiques
[Guo ef all92l Y. Q. Guo, J. L.Batoz, M. Elmouatassimet J. M. Detraux, On the estimation

of thickness strain in thin car panels by the inverse approach, NUMIFORM'92, Chenot,
Wood, Zienkiewicz (Eds),Balkema,Rotterdam,1992.
tHadji/951 S. Hadji, Méthode de résolution pour les fluides incompressibles.Thèse de

doctoratde I'Université de Technologiede Compiègne,1995.
[Halphen et Salençon/87] B. Halphen et J. Salençon,Elasto-plasticité,Pressesde I'Ecole

NationaledesPontset Chaussées,Paris, 1987.
[Hacquin/96] A. Hacquin, Modélisation thérmomécaniquetridimensionnelledu laminage,

Thèsede doctoratde I'Ecole Nationale SupérieuredesMines de Paris, 1996.
tHilVs0l R. Hill. The mathematicaltheory of plasticity, ClarendonPress,Oxford (G. B.),

1950.
[Hunek/93] I. Hunek, On a penality formulation for contact-impactproblem,s.Computersand

Structures,48, pp.: I 93-203,1993.
[Imazatene/01] A. Imazatene, Méthodes numériques pour les problèmes non linéaires :

Application à l'élastoplasticité en transformationsfinis, Thèse de doctorat de I'université de
Metz. En cours.
èmeédition,
[Imbert/9l1 J. F. Imbert, Analyse des structures par élémentsfinis, 3
CEPADUES, Toulouse,1991.
[Jean et Touzot/88] M. Jean et G. Touzot, Implementationof unilateral contact and dry

friction in computer codes dealing with large deformation problems, Journal de Mécanique
Théoriqueet Appliquée,Specialissue,I,7,pp:.145-160,1988.
[Jean/93f M. Jean, Simulation numérique des problèmes de contact avec frottement,

Matériaux et Techniques,l-2-3, pp.: 22-32, 1993.
[Ju et Taylor/881 J. W. Ju, et R. L Taylor, A perturbed lagrangian formulation for thefinite

elementsolution of non linearfrictional contactproblem,s,Journal de Mécanique Théorique et
Appliquée,Specialissue,l, 7, pp.:1-14,1988.
[Kallala/93] C. Kallala, Approche géométrique du contact entre solides par la méthode des
élémentsfinis dans un cadre grandes déformationsélastoplastiques,Thèse de doctorat de
L'ENSAM, Paris, 1993.
[Kikuchi et Oden/88] N. Kikuchi et J. Oden, Contqct Problems in Elasticity: A Study of

Variational Inequalities and Finite Element Methods, SIAM StudiesIn Applied Mathematics,
1988.
[Klarbringllsfil A. Klarbring, A mathematicalprogramming to three dimensional contact
problem with friction, Computer Methods in Applied Mechanicsand Engineering,58, pp.:
175-200,1986.
-tr4-
R4férences bi blio qraphiques
[Kobayashi et Kim/78] S. Kobayashi et J. Kim, Deformation analysis of axisymmetric sheet

metalforming processesby the rigid plastic finite elementmethod, Mechanics of SheetMetal
Forming, (Eds) D. P. Koistinen and N. M. Wang, Plenumpress,New York, pp;341,1978.
[Koiter/45| W. T. Koiter, On the stability of elastic equilibrium. Thèse de doctorat de

I'Université de Delft, English translation:NASA TechniquesTranslations,F. 10,883, 1967,
t945.
dans la mise en
[Lachat et Ravier/?qt l. C Lachat et P. Ravier, Rechercheet développement
æuvre des aciers,4'-" Congrèsde Calcul des Structures,D. Guédra-Degeorges, P. Ladevèze
et M. Raous,18-21mai, Giens,pp.:4l-49,1999.
[Ladevèze/8s] P. Ladevèze, Sur une famille d'algorithmes en mécanique des structures,

ComptesRendusAcad. Sc.Paris,300,2,pp.:41-44,1985.
[Ladevèze/g| P. Ladevèze, New advances in the large time increment method, New
advancesin Computational Structural Mechanics,P. Ladevèze - O.C. Zienkiewicz (Eds),
Elsevier,pp.:3-21,1991.
[Ladevèzel96l P. Ladevèze,Mécanique non linéaire des structures: Nouvelle approche et

méthodesde calcul non incrémentales,Hermès,Paris, 1996.
[Lee et KobayashilT3l C. H. Lee et S. Kobayashi, New solution to rigid-plastic deformation

problems using a matrix method, Journal of EngineeringIndustries,ASME, 95, pp.: 865,
1973.
[Lebon/881 F. Lebon, Résolution numérique des problèmes de frottement de Coulomb.

Accélération de la convergencepar une méthodemulti-grilles internes, Thèse de doctorat de
I'université de ProvenceAix-Marseille, 1988.
[Lemaitre et Chaboche/85]J. Lemaitre et J. L. Chaboche,Mécaniquedes matériauxsolides,

Dunod,Paris,19985
[Lochegnies et allï9l D. Lochegnies, J. Oudin, J. Picard, J.C Jelin, J.M. Rigaut et Y.

Ravalard, Frictional contact modeling in viscoplasticity and large deformation by finite
element method, Fifth lnternational Symptium on Numerical Method in Engineering,
Lausanne,1989.
[Nagtaat et Taylor/gl] J. C. Nagtaal et L. M. Taylor, Comparisonof implicit and explicit

finite elementmethodsfor analysis of sheetforming problems, VDI BERICHTE, 894, pp.:
705-725,t991.
[Najah et allgïl A. Najah, B. Cochelin, N. Damil et M. Potier-Ferry,A critical review of

asymptotic numerical methods, Archives of Computational Methods in Engineering, 5, pp.:
31-50,1998.
[Noor et Peters/80] A. K. Noor et J. M. Peters,Reduced basis techniquefor nonlinear

analysisof structures,AIAA Journal, 18, 4, Article No 79-0747R,1980.
1 1 5-
Ré-fërences bib li o qraphiques
[Noori8l] A. K. Noor, Recent advances in reduction methodsfor nonlinear problems,

Computersand Structures,13,pp.: 31-44,1981.
[Nour et Peters/8l] A. K. Noor et J. M. Peters,Tracingpost-limit-pointpaths with reduced

basistechnique,ComputerMethodsin Applied Mechanicsand Engineerin5,2S,pp.:217-240,
198r.
[Oden et Pires/83] J. T. Oden et E.B. Pires, Non local and non-linearfriction laws and
variational principles of contact problems in elasticity, Joumal of Applied mechanics, 50,
pp.:67-75,1983.
[Oden et Pires/88] J. T. Oden et E.B. Pires,Analysis of contactproblems with friction under

oscillating loads, Computer Methods in Applied Mechanics and Engineering, 39, pp.:337-
362,1988.
[Oden/88f J. T. Oden, New models of friction for non linear elastodynamicsproblems,

Journalde MécaniqueThéoriqueet Appliquée,Specialissue,l ,7,pp.;47-54, 1988.
[Onâte et Saracibar/921E. Onàte et C. Agelet de Saraclbar,Alternativesfor finite element

analysis of sheet metal forming processes, Proceedingsof Numerical Method in Industrial
Forming Processes, Chenot,'Wood,Zienkiewicz (Eds),pp.: 78-88,1992.
[Papadopoulos et Taylor/93] P. Papadopouloset R. L. Taylor, A simple algorithm for

threedimensionalfinite elementanalysis of contactproblems, Computersand Structures,46,
p p . :l l 0 7 - 1 1 1 8 , 1 9 9 3 .
[Peeters et Konter/gz] F. J. H. Peeters et A.W.A. Konter (MARC Analysis Research

Corporation-Europe), Numerical simulation of industrial metal forming processes,
Numiformeg2; Sophia-Antipolis,I 4- 18 September1992.
[Picart et all99] P. Picart, O. Ghouati,N. Boudeau,J. C. Gelin, Optimisationdesparamètres

de procédés de miseforme des matériaux avec minimisation de I'endommagementet de la
striction localisée,4"-" Congrèsde Calcul des Structures,D. Guédra-Degeorges, P. Ladevèze
et M. Raous,l8-21 mai, Giens, pp.: 355-360,1999.
[Pinto/90f Y. Pinto, Contact et frottement en grande déformation plastique, Application au

serre-flan en emboutissage,Thèse de doctorat de I'université de ProvenceAix-Marseille I,
1990.
[Potier-Ferry/87] M. Potier-Ferry,Foundations of elastic post-buckling theory, Buckling and

Post-buckling,LectureNotes in Physics,288, pp.:1-82,Springer-Verlag,Berlin, 1987.
[Potier-Ferry et all97l M. Potier-Ferry,N. Damil, B. Braikat, J. Descamps,J. M.Cadou, H.

L. Cao et A. Elhage-Hussein, Traitement des fortes non-linéarités par la méthode
asymptotique numérique,C. R. Acad. Sci. Paris,324, SérieII b, pp.: l7l-177 , 1997.
[Rachik et Roelandt/gSl M. Rachik et J. M. Roelandt,Schémaexplicite d'ordre élevépour

les problème non linéaires de calcul de structures, Revue Européennedes Eléments Finis, 7,
4 pp.: 401-420,1998.
- l16 -
Référencesbiblioetsplhue;
[Ramm/8l] E. Ramm, Strategiesfor tracing the nonlinear response near limit points. Non-
'Wunderlich,
linearfinite elementanalysis in structural mechanics,(Eds) E. E. Stein and K.J.
Bathe,Springer-Verlag,Berlin, 1981.
[Raous et al l88l M. Raous,P. Chabrandet F. Lebon,Numerical methodforfrictional contact

problem and applications,Journalde MécaniqueThéoriqueet Appliquée,Specialissue, 1, 7,
p p . :1 l l - 1 2 8 , 1 9 8 8 .
[Riks/S4f E. Riks, Some computational aspects of the stability analysis of nonlinear

structures,ComputerMethods in Applied Mechanicsand Engineering,47,pp.:219-259,1984.
[Signorini/30] Signorirri, Sulle deformazioni termoelastichefinite, 3'd International Congress

of Applied Mechanics,2,pp.:80-89,1930.
[Thompson et Walker/68] J. M. T. Thompson and A. C. Walker, The non-linear

perturbation analysis of discrete structural systems, lnternational Joumal of Solids and
Structures,4, pp.:757-768, 1968.
[Simo et RifaVg0] J. C. Simo et M. S. Rifaï, A class of mixed assumedstrain methodsand

method of incompatible modes, lnternational Journal for Numerical Methods in Engineering,
37, pp.;I 595-1636,1990.
[Simo et Laursenlgzl J. C. Simo, T. A. Laursen,An augmentedLagrangian treatment of

contactproblems involvingfriction, Computersand Structures,1, I ,pp.:97-1t6,1992.
[Tang/8lf S. C Tang , Large elastoplasticstrain analysis of flanged holeforming, Computers

andStructures,13,pp.:363,1981.
[Tifler/83] W. Tiller, Rational B-Splines for curve and surfoce representation, IEEE,
ComputerGraphicsand Applications,1983.
[Tri/96] A. TRI, Méthode asymptotiquenumériql.epour lesfluides visqueuxincompressibles

et la détection de la bifurcation de Hoff,Thèse de doctorat de I'universitéde Metz,1996.
lTri et all96l A. Tri, B. Cochelin et M. Potier-Ferry,Résolutiondes équationsde Navier

Stokeset détection des bifurcations stationnairespar une Méthode AsymptotiqueNumérique,
RevueEuropéennedesElémentsFinis, 5, 4, pp.:415-442,1996.
[Vannucci et all9ïl P. Vannucci, B. Cochelin,N. Damil et M. Potier-Ferry,An asymptotic

numerical method to compute bifurcating branches, International Journal of Numerical
Methodsin Engineeing, 41,pp.: 1365- I 3 89, I 998.
'Walker,
[Walker/69] A. C. A non-linearfea of shallow circular arches,lnternationalJournal
of Solids and Structures,5, pp.:97-102,1969.
[Wang/841 N. M. W*g, A rigid-plastic rate-sensitivefinite elementmethodfor modelling

sheet metal forming processes,Numerical analysis of Forming Processes,(Eds) J. F. T.
Pitman,O. C. Zienkiewicz, R. D. Wood et J. M. Alexander,John Wiley & SonsLtd, London,
p p . :l l 7 - 1 6 4 , 1 9 8 4 .
-tt7 -
R4férences bib lio grap hiques
[Wang et Budiansky/78] N. M. Wang et B. Budiansky,Analysis of sheetmetalforming by a

finite elementmethod,Joumalof Applied Mechanics,ASME, 45, pp.: 73,1978.
['Wagoner a all88l R. H. Wagoner, E. Nakamachi et Y. Germain, Experience with explicit

and implicit finite elementprogramsfor sheetforming analysis,Modelling of Metal Forming
Processes,(Eds) J. L. Chenotand E. Onate,Kluwer Academicpublishers,pp.: 131, 1988.
tWifVT6l A. S. Wifi, An incremental completesolution of the stretch forming and deep

drawing of a circular blank using an hemispherical punch, lnternational Journal of
MechanicalSciences,18,pp.: 23,1976.
[Wood et all89l R. D. Wood, J. Bonnet et A. H . S. Wagadipura,Simulation of the

superplastic forming of thin sheet components using the finite element method,
NUMIFORM'89, (Eds) Thompson,Balkema,Rotterdam,1989.
lZthrounilgsl H. Zaltroln| Méthode AsymptotiqueNumérique pour les coques en grandes

rotations,Thèsede doctoratde I'universitéde Met2,1998.
[Zahrouni a atDSl H. Zahrouni, M. Potier-Ferry,H. Elasmar et N. Damil, Asymptotic

numerical methodfor non-linear constitutive/aws, Revue Européennedes ElémentsFinis, 7,
pp.: 841-869,1998.
[Zahrouni et atl99'lH. Zahrouni, B. Cochelin et M. Potier-Ferry,Computingfinite rotations

of shells by an asymptotic numerical method, Computer Methods in Applied Mechanics and
Engineering, 175,pp.: 7l-85, 1999.
lZhang et atl86l W. C. Zhang, R. D. Wood et O. C. Zienkiewicz, Superplasticforming

analysis using a finite element viscousflow formulation, Aluminium Technology'86, The
instituteof metal,London,4,pp.:111.4-111.6, 1986.
[Zienkiewicz et atl78] O. C. Zienkiewicz, P. C. Jean et, E. Onate, Flow of solids during

14, pp.:
forming and extrusion,some aspect of numerical solutions, Solids and Structures,
15.1978.
[Zienkiewicz et Taylor/9l1 O. C. Zienkiewicz et R. L. Taylor, The Finite ElementMethod,

Mc Graw-Hill, 4th Edition, London- New-York, 1991.
lZhong!93| Z.H. Zhong, Finite elementprocedure for contact-impactproblems. Oxford

University press,Oxford, 1993.
-118-
Simulation de I'emboutissageà froid par une Méthode Asymptotique Numérique
L'emboutissageest un procédé de mise en forme des tôles minces, qui permet de réaliserde
nombreux objets ou pièces, en particulier pour I'automobile. La simulatiou complète de pièces
complexesrestedifficile à réaliseren raison de la ldngueurdes calculsà effectuer,qui est due à la taille
des problèmes à résoudre (nombre de degrés dJ tberté) et au couplage de multiples non-linéarités
(grandes transformations,comportementplastique,contact, frottement).L'objectif de cette thèse est
d'alléger ces calculs en utilisant une nouvelle *éthod" de résolutionappeléeMéthode Asynrptotique
Numérique. Cette méthode consiste à transformer les problèmes non-linéairesà résoudre en utre
successionde problèmeslinéairesqui admettenttous le même opérateurtairgent.Dans cette thèseon a
adoptéune formulationde coquebien adaptéeau phénomènesde grandestransformatiou.Att niveau de
la loi de comportement,la déchargeélastiquen'estpas prise en compte,ot'.a adoptédonc la théorie de
la plasticité en déformationtotale. Pour réaliser des calctrlsde sériesentières,il a été nécessairedc
régulariserla loi de comportelnentainsi quc la loi de contactunilatéral. Le couplagedes dilferentes
linéarité n'a pas altéré I'ef{lcacitéde la IVIAN. Les résultatsobtcuus monffcnt bien la précision des
Cerésolutionclassiques.
algorithmeset la rapiditéde la métlrodeutiliséeen la comparantaux métho<les
Mots clé
- plasticitj - coutîct - méthodcde perfurbation- uÉthode clesélémentsfinis.
Emboutissage
Sirnulationof thc cokl dccp drarving by an As1'rnptoticNumcrical Nlethocl
'Ihe
dcep rlrau'ing is a fonning l)t'ocesslargely used in the inclustry,especially in the car'
manufechrres.Thc nruuelicrl simulation o[ this kind of ploccssesncctls tttttclt CPU tinre bcc:tuseof
anclthc constihrtivelaw.'['he ain-r
du: to the !;eonrtr-v,contactphenomenon
the severalnonlincarities
reù,:'.,.ethc sonrputing tinre lrsitr,: ltn llt.rrnative nrctltcil cir-llecl"-!.,:',trrptotic
of this thesi; is 1,.1,
Numer;cal Vlethod". This lattei allorvsus to searcirsolltiou br':nchesiIr t]tc t-onittfçç,n x lcries aucl
I tnns:brms the nonlinerr problenu into a successionoî rvell posedlinear problernsinvolving ths same
tangent operatol. In this tiresis a shell formulatiou well adopted to large transfonnation phenomenonrs
adopted. As the elastic unload is considered,the theory of plasticity rvith the total deforntationis
adopted. To apply a perturbationtechnique,a regularisationof the constirutivelarv as rvell as the
ccnrscrconditionshas been used.Comparedto the classicaliterativemethod,the A.N.M. needsahvays
lcss mctrix decompositions anci less computing time. The reason of this peifotmance is the abilitv of
the method to deconrposeonly one stiffness matrix by step and to adjust automatically the step size of
the local nonlinearify encountered.For the class of problems presentedin this thesis, the ANM is
efficient, reliable and easyto perform.
Kev rvords
Deep drawing - plasticity - contact - Perturbation techniques- finite element method.

Abichou Hammadi SMZ0110

Transféré par

Droits d'auteur :

Formats disponibles

Abichou Hammadi SMZ0110

Transféré par

Informations du document

Description originale:

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Abichou Hammadi SMZ0110

Transféré par

Droits d'auteur :

Formats disponibles

Simulation de l’emboutissage à froid par une méthode

To cite this version:

HAL Id: tel-01775454

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

Ce document est le fruit d'un long travail approuvé par le jury de

Il est soumis à la propriété intellectuelle de l'auteur. Ceci

D'autre part, toute contrefaçon, plagiat, reproduction illicite

Code de la Propriété Intellectuelle. articles L 122. 4

INSTITUT SUPERIEURDE GENIE MECANIQUE

Pour I'obtention du gradede :

J. C. Gelin LMARC & ENSMM Besançon Rapponeur

ISGMP',Ile du saulry' 57M5METZ Cedex01

Je tiens égalementà remercier MessieursM. Gueury,Y. Chevolier, H. Ben

le remerciemesomieset omisdu L.P.M.M.pour leur oide el leur soutien

,4 ^o fn**n, mzr parenù et a'hnh rno fo*i//n, en Junirin û nn 7ron n...

CHAPITRE I Etude biblioeranhique 9

1.2. Intérêt d'une simulation numérique du procédé T2

1.3. Modélisations mécaniquesdu procédé de I'emboutissage l5

1.3.1.Modèle rigide Plastique 15

1.4. Eléments finis utilisés dans les simulations numériques de I'emboutissage t7

l.4.2.Approches bidimensionnelleset approchessolides l8

1.5. Problème de contact et de frottement l9

1.5.1.Positiondu problème de contact 20

f .5.3. Lois de frottements utilisées dans les simulations numériques 22

1.5.3.c.Modèle defrottement non local 24

1.5.4.Méthodes de traitement des problèmes de contact et de frottement 25

1.5.4.b.Méthode desMultiplicateurs de Lagrange 27

1.5.4.d.Méthode de Lagrange augmentée 29

1.5.4.e.Méthode des élémentsd'interfoce 30

1.6. Méthodes de résolution numériques 31

1.6.2.Méthodes Asymptotiques Numériques 32

2.1. Formulation du problème 44

2.1.4.Etude du problèmede contact 48

2.1.4.a. Conditions de contact 49

2.1.4.b.Régularisation de la loi de contact 50

2.1.4.c. Géométriede contact en 3D 5l

2.2.2.Les approximantsde Padé 56

CHAPITIIE 3 Application à la simulation de I'emboutissaeehémisphérioue 60

3.1. Simulation de I'emboutissagehémisphérique 61

3.1.2.Traitement des problèmes non-linéaires dans Abaqus 64

3.2. Etude de loordre optimal de troncature des séries 69

3.3. Etude de I'influence des paramètres de régularisation 7l

3.3.1.Influencede la régularisationde la loi de comportement 7l

3.3.2.Influencede la régularisationde la loi de contact 72

3.4. Efficacité de I'algorithme 75

3.5. Etude de I'amincissementde la tôle 77

3.6. Etude de la qualité de la solution 79

CHAPITRE 4 Adoucissementde la loi de contact 83

CHAPITRE 5 Généralisation des surfacesdes outils 95

5.1. Méthodes utiliséespour la représentation de surfacesde contact 96

5.2. Représentationdes surfacespar les fractions rationnelles 99

Conclusions généraleset perspectives r07

La réalisation par emboutissage, des pièces à parois minces, est couramment

La résolution numérique d'un problème de mise en forme consiste à déterminer les

Les logiciels actuels de simulation numérique utilisent le plus souvent des

Ces dernières années, plusieurs travaux ont été consacrésà I'amélioration de la

Dans le premier chapitre, on présenteune étude bibliographique sur I'emboutissage

Dans le deuxième chapitre, on présente une modélisation du problème de contact