These E ROUX C
These E ROUX C
These E ROUX C
Doctorat ParisTech
THÈSE
pour obtenir le grade de docteur délivré par
Emile ROUX
le 07 juillet 2011
Assemblage Mécanique :
Stratégies d’optimisation des procédés et d’identification des
comportements mécaniques des matériaux
Jury
M. Jacques BLUM, Professeur, Laboratoire J.A. Dieudonné, Université de Nice Sophia Antipolis Président T
M. Michel BRUNET, Professeur, LamCoS, INSA Lyon Rapporteur
M. Fabrice SCHMIDT, Professeur, ICA-Albi, Ecoles des Mines d'Albi Rapporteur
H
M. Mohamed RACHIK, Maître de Conférences - HDR, Laboratoire Roberval, UTC Compiègne Examinateur
È
M. Julien MALRIEU, Ingénieur, CETIM Invité
M. François BAY, Professeur, CEMEF, Mines ParisTech Directeur de thèse S
M. Pierre-Olivier BOUCHARD, Maître de Conférences - HDR, CEMEF, Mines ParisTech Co-directeur de thèse
E
MINES ParisTech
Centre de Mise en Forme des Matériaux
Rue Claude Daunesse B.P. 207, 06904 Sophia Antipolis Cedex, France
2
Remerciements
Je tiens ensuite à remercier Fabrice Schmidt et Michel Brunet qui m'ont fait l'honneur de
rapporter mon travail de thèse. J'adresse également mes remerciements à Jacques Blum,
Mohamed Rachik et Julien Malrieu qui ont examiné mon travail.
Mes remerciements vont ensuite vers toutes les personnes que j'ai pu rencontrer au
laboratoire et au Cetim, personnes qui ont rendu ces trois années très enrichissantes grâce
aux collaborations et discussions que nous avons pu avoir. Notamment, mes collègues de
projet Rachid El Kahoulani, Mathieu Gaurin et Julien Malrieu. Merci aussi à toutes celles et
ceux qui m'ont permis d’égayer les soirées, les Week End et les pauses café durant ces
années. Je remercie plus particulièrement Sabine, Franck, Thomas, Nicolas, Lionel, Aurélien,
Benoit, Hicham, Nadine, Benjamin, Sylvie, Ziad, ...
Je souhaite aussi exprimer toute mon amitié à des personnes que j'ai rencontré au Cemef
et qui sont devenues des Amis: Raph, Pam, Greg, Larbiao et Tof. Je les remercie aussi bien
pour nos séances de travail que de délire, sans quoi ces années auraient été bien ternes,
Merci à vous!
Je profite aussi de ces quelques lignes pour remercier une personne que j'estime
beaucoup, qui se reconnaitra sûrement, qui m'a donné l'envie et la confiance en moi
nécessaire pour me lancer dans cette aventure, Merci!
Enfin je remercie ma famille de m'avoir soutenu et toujours encouragé durant toutes ces
années d'études.
3
Sommaire
Introduction................................................................................................................................. 11
2 Bibliographie.................................................................................................................................. 18
2.1 Généralités sur l’optimisation.................................................................................................................18
2.2 Analyse de sensibilité par plans d’expériences...............................................................................28
2.3 Méthodes d’optimisation pour fonctions objectif coûteuses ...................................................34
2.4 Bilan....................................................................................................................................................................39
5 Conclusion ....................................................................................................................................... 55
5
1.4 Architecture logicielle retenue – environnement objet...............................................................61
1.5 Gestion des calculs – exploitation des ressources de calcul ......................................................62
1.6 Evolutivité de la plateforme.....................................................................................................................63
4 Conclusion ....................................................................................................................................... 89
6
4 Conclusion .................................................................................................................................... 170
Références ..................................................................................................................................209
7
Table des abréviations et des sigles
AE Algorithmes évolutionnaires
AG Algorithmes génétiques
BFGS Méthode de minimisation Broyden-Fletcher-Goldfarb-Shanno
CL Constant Liar
CV Cross Validation
EF Eléments finis
EGO Efficient Global Optimisation
EGO-p Algorithme EGO parallèle
EI Expected Improvement
KB Kriging believer
LB Lower confidence Bound
MDC Mesures de champs
MLE Maximum Likelihood Estimation
MOOPI MOdular software dedicated to Optimisation and Parameters Identification
PE Plans d'expériences
SE Stratégies d’évolution
SE-Meta Stratégies d’évolution assistées par Métamodèle
9
Introduction
Introduction
L'un des grands défis de l'industrie est d'optimiser l'ensemble du cycle de vie d'un produit,
depuis la conception jusqu’à la tenue mécanique en service. L'objectif est de réduire les coûts
au sens large du terme, c'est-à-dire aussi bien le coût financier que le coût environnemental.
La phase de conception est l'une des étapes clés pour atteindre ces objectifs. L'utilisation
d'outils de conception assistés par ordinateur est une réponse très efficace à ces
problématiques.
Pour répondre à ces exigences de plus en plus élevées dans l'industrie (aéronautique,
automobile, de l'électroménager, du transport ferroviaire…), les concepteurs ont recours à
l'utilisation de matériaux de natures différentes et sont confrontés aux problématiques
d'assemblage de ces matériaux. Les techniques d’assemblage par soudage ou par collage sont
souvent retenues, mais ne permettent pas de répondre à toutes les problématiques rencontrées
dans l’industrie. Différentes techniques d’assemblage par déformations plastiques peuvent
alors être envisagées. Ces procédés d'assemblage par déformations plastiques regroupent des
procédés très anciens tels que le rivetage mais aussi des procédés plus innovants tels que le
clinchage, le sertissage, ou encore le rivetage auto-poinçonneur… Ces procédés présentent
l’avantage de pouvoir créer un assemblage entre des matériaux ayant de mauvaises propriétés
de soudabilité ou de collage. La plupart de ces procédés présentent de bonnes propriétés de
tenue mécanique, principalement en cisaillement, et sont fiables dans le temps. Certains
d’entre eux, tel que le clinchage, sont également intéressants d’un point de vue financier,
puisqu’ils ne nécessitent pas d’apport de matière supplémentaire.
La maîtrise et la bonne connaissance de ces procédés constituent donc des points
primordiaux pour pouvoir les intégrer de manière efficace et optimale dans les structures
mécaniques industrielles. Le développement d'outils numériques capables de modéliser les
assemblages par déformation plastique est par conséquent un point clé pour l'aide à la
conception et à l'innovation. Ce type d'outils vient compléter les outils et le savoir déjà
existant, basé sur des observations empiriques et l'expérience des concepteurs.
C'est dans ce contexte d'aide à la conception et à l'innovation que le projet MONA LISA
(Modélisation et Optimisation Numérique des Assemblages, Logiciel Intégré de Simulation
des Assemblages) a été construit. Il implique 3 partenaires : le CETIM (Centre technique des
industries mécaniques), le CEMEF (Centre de mise en forme des matériaux) et Transvalor.
Les objectifs du projet peuvent se décliner en deux volets. Le premier est de développer les
méthodes numériques nécessaires à la modélisation et l'optimisation des technologies
d’assemblage par déformation pastique. Le second est de proposer un ensemble de logiciels
11
d'aide à la conception et à l'optimisation des assemblages, intégrant l'ensemble de ces
développements numériques.
Problématiques de la thèse
Cette thèse est l’une des quatre thèses effectuées dans le cadre du projet MONA LISA. La
problématique traitée ici est double :
- Mise en place d’outil numériques permettant la caractérisation des propriétés mécaniques
des matériaux mis en jeu dans l'assemblage.
- Développement d’un outil d'aide à la conception et plus particulièrement à l'optimisation
des procédés d'assemblage par déformation plastique pour améliorer leur tenue en service.
L'ensemble des travaux réalisés au sein du projet MONA LISA permettra, à terme, de
mettre en place une chaîne de modélisations complète depuis la mise en place de l'assemblage
jusqu’à la tenue en service. Ces modélisations permettront donc d'analyser une configuration
d'assemblage donnée, et d’évaluer l'influence des paramètres du procédé sur la tenue
mécanique en service. Mais pour exploiter pleinement la puissance d'un tel outil il est
incontournable d'utiliser des méthodes d'analyse et d'optimisation, permettant de proposer de
nouvelles solutions d'assemblage, et d'optimiser un assemblage par rapport à une sollicitation
mécanique définie. La première problématique posée est donc de proposer une solution
d'optimisation et d'analyse des assemblages par déformation plastique exploitant les outils de
modélisation développés dans le cadre du projet. En filigrane de cette problématique apparait
une problématique connexe liée au temps de calcul. En effet les calculs permettant la
modélisation de la mise en place et de la tenue mécanique d'un assemblage par déformation
plastique sont relativement longs (plusieurs heures). Si l'on souhaite tester, de manière
automatique, plusieurs configurations d'un assemblage pour aboutir à une solution optimale
en un temps raisonnable, il est nécessaire d'utiliser des méthodes intégrant cette
problématique de temps.
Outre les méthodes numériques utilisées, la précision des résultats est liée à la bonne
connaissance du comportement des matériaux composant l'assemblage. La modélisation des
procédés d'assemblage par déformation plastique met en jeu de grandes déformations
plastiques ainsi que des processus d'endommagement ductile dans les différentes pièces
composant l'assemblage. On se posera donc comme seconde problématique l'identification de
ces paramètres matériaux et notamment l'identification des paramètres matériaux liés aux lois
d'évolution de l'endommagement ductile. La prédiction fiable de l'endommagement étant un
point primordial pour bien évaluer la tenue en service d'un assemblage.
12
Introduction
analyse inverse et à l'optimisation. En plus de la réflexion sur la mise en place d'une structure
logicielle pérenne, ce chapitre fait une large place à la validation et à l'évaluation de
l'efficacité de la plateforme et de l'algorithme d'optimisation proposé.
Les trois chapitres suivants exploitent cette plateforme.
Le chapitre 3 propose une procédure d'optimisation des assemblages par déformation
plastique. En préambule de ce chapitre, une brève description du modèle de comportement du
matériau est faite : le modèle retenu pour la modélisation des assemblages par déformation
plastique est constitué d’une loi de comportement élasto-plastique couplée à un modèle
d’endommagement ductile de Lemaitre. La procédure mise en place est illustrée par
l'optimisation du procédé de clinchage. Cette procédure se décompose en trois phases : la
validation de la modélisation du procédé, l'identification des paramètres procédé les plus
influents, et l'optimisation du procédé.
Le chapitre 4 propose une méthodologie d'identification par analyse inverse des
paramètres matériaux, analyse inverse basée sur l'exploitation de la plateforme MOOPI. Ce
chapitre est dédié à l'identification des différents paramètres du modèle de comportement
retenu sur la base d'essais de traction classique.
Puis, dans le chapitre 5, une méthodologie pour intégrer les mesures de champs dans le
processus d'identification est mise en place. L'objectif de ce chapitre est de montrer l'intérêt
d'exploiter des observables plus riches pour identifier les paramètres matériaux de manière
plus pertinente.
En conclusion de ce manuscrit, le bilan des travaux est dressé et une discussion sur les
ouvertures et perspectives envisagées est proposée.
13
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
Chapitre 1
Optimisation dans le cadre de la modélisation par éléments
finis des procédés d'assemblage
15
Chapitre 1
1.1 Définition
Les méthodes d'assemblage par déformation plastique regroupent les méthodes
d'assemblage, à froid, de produits minces dont la tenue mécanique est assurée par la
déformation et l'interpénétration des composants de l'assemblage. Ces assemblages peuvent
être créés sans apport de matière (clinchage Figure 1-1 c et d) ou avec apport d'un composant
supplémentaire (rivetage Figure 1-1 b, rivetage auto-poinçonneur Figure 1-1 b)
(a) (b)
(c) (d)
Figure 1-1 : Différentes techniques d'assemblage par déformation plastique : (a) rivetage auto-
poinçonneur (coupe), (b) rivetage, (c) clinchage à point rond, (d) clinchage à point carré
16
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
Figure 1-2 : Résultats de l'optimisation d'un assemblage par rivet auto-poinçonneur [Fay08]
Les résultats obtenus montrent l'intérêt d'optimiser le procédé d'assemblage. Sur cet
exemple cependant, une seule variable a été considérée et l’optimisation a été réalisée de
manière manuelle. Une optimisation automatique de la chaine virtuelle de simulation (pose du
point d’assemblage et modélisation de la tenue mécanique) constituerait un apport majeur en
terme d’efficacité et d’utilisation dans un cadre industriel. L'objectif final serait donc de
réaliser une boucle d'optimisation automatique (Figure 1-3) capable de travailler sur plusieurs
paramètres du procédé afin d'améliorer au maximum la tenue mécanique.
Simulation Simulation
Pose de l’assemblage Tenue mécanique en service
Optimisation
Dans la suite de ce chapitre une revue des différentes méthodes pour optimiser et
améliorer la tenue mécanique des assemblages sur la base d'expériences numériques est
présentée.
17
Chapitre 1
2 Bibliographie
Dans la plupart des cas, les paramètres d'entrée d'un processus sont très nombreux. Pour
un point d'assemblage de deux tôles par rivetage, on peut lister : les diamètres des perçages
des deux tôles, le diamètre du rivet, la longueur du rivet, la forme de la tête du rivet, la vitesse
de la bouterolle, l'effort de pose, la forme des outils... Tous ces paramètres ont une influence
plus ou moins importante sur la tenue mécanique de l'assemblage et peuvent donc être
considérés comme des variables d'optimisation.
d) Domaine d'optimisation
La définition des variables d'optimisation pousse naturellement à définir le domaine
d'optimisation, c’est-à-dire l'espace de variation des paramètres x . On le notera X . Dans le
cadre de ce manuscrit on travaillera exclusivement sur des domaines d'optimisation continus.
x∈ X
(1.3)
X ⊂ ℜn
Il est à noter qu'il existe des méthodes adaptées aux domaines non continus, comme par
exemple les plans d'expériences (présentés au paragraphe 2.2) qui peuvent prendre en
considération ce type de variables.
Dans le cadre de l'optimisation des assemblages, la définition des bornes de variation des
paramètres d'optimisation permet de définir ce domaine d'optimisation.
18
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
e) Contraintes
Les contraintes peuvent être classées en deux catégories : les contraintes explicites et les
contraintes implicites.
Les contraintes explicites sont directement appliquées aux paramètres d'optimisation ; le
jeu de paramètres x peut être déclaré non conforme sans même évaluer la solution complète.
Par exemple pour un assemblage par rivetage le rivet devra avoir un diamètre inférieur ou
égal au diamètre de perçage.
Les contraintes implicites ne sont, quant à elles, appliquées qu'après évaluation complète
d'une solution.
Ces contraintes peuvent être représentées mathématiquement sous forme d'égalité ou
d'inégalité entre les variables d'optimisation.
Algorithme Evaluation de
d'optimisation la solution
f ( x), g ( x), h ( x)
19
Chapitre 1
a) Principes généraux
Supposons que la fonction coût soit continue et dérivable sur l’ensemble de l’espace de
recherche, et notons ∇f ( x ) le gradient de la fonction coût en x .
La condition nécessaire d’optimalité s’écrit :
∇f ( x ) = 0 (1.5)
Si f ( x ) est deux fois différentiable il est possible d'écrire la condition suffisante
d'optimalité :
f ( x) = 0
∇
2 (1.6)
∇ f ( x ) définie positive
Les méthodes à direction de descente ont pour objectif de calculer le vecteur x satisfaisant
la condition nécessaire (1.5).
k +1 (1.9)
x = x + λ d , ∀k ≥ 0
k k k
Cette présentation générale fait apparaître les deux paramètres de base des algorithmes à
direction de descente : la direction de descente d et le pas de descente λ .
20
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
• La dérivabilité de la fonction coût par rapport aux variables d'optimisation doit être
vérifiée, la fonction ne doit pas être bruitée. Dans le cadre de calculs éléments
finis, cela peut poser problème lorsque le calcul de la fonction coût est sensible à la
finesse du maillage ou au remaillage.
• La solution obtenue vérifie la condition nécessaire d’optimalité (1.5), cette
condition permet uniquement de conclure que la solution proposée est un
minimum local et non global sur le domaine. Cette solution dépend fortement du
point de départ x 0 (1.9).
21
Chapitre 1
( ) ( )(
Q ( x ) = f x k + ∇f x k . x − x k + ) 1
(
t
) ( )(
x − x k .∇ 2 f x k . x − x k ) ∇2 f ( x ) f
de
2
On choisit x k +1 qui minimise Q( x ) , on obtient alors la direction de
descente grâce à la condition d’optimalité.
Gauss- Extension de la méthode de Newton basée sur une approximation du Largement
Newton /
hessien de
f dans le cas où la fonction coût s'exprime sous forme de exploitée pour
Levenberg- l'identification
moindres carrés [Min98].
Marquardt de paramètres
( ) ( )
−1
Approximation itérative du hessien de f .
Quasi- Efficacité
Newton
d k =- Hɶ k . ∇f x k
Les méthodes les plus populaires pour construire cette approximation
Hɶ k sont DFP ([Fle63]), et BFGS [Bro70].
22
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
40
30
20
10
f
x’
0
0 x
0 d’
d x
l
−10
−20 x
g
−30
0 2 4 6
x
Pour dépasser ces problèmes de minima locaux et de dérivabilité, il existe les algorithmes
dits globaux. Ils font l'objet des paragraphes suivants.
23
Chapitre 1
24
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
Population initiale
Évaluation de la population
Population parents
Opérateur génétique
(sélection, croisement, mutation)
Actualisation de la population
non
Test de
convergence
oui
Sortie : extraction du meilleur individu
- Initialisation de la population
La population initiale est générée dans le domaine d'optimisation. Le nombre d'individus
est fixé par l'utilisateur. Cette génération peut être soit aléatoire, soit suivre des règles
d'échantillonnage comme les méthodes Latin Hyper Cube [Tan93].
- Évaluation de la population
L’évaluation de la population est une étape qui peut devenir coûteuse si le temps
d'évaluation de la fonction coût est important.
- La sélection
La sélection permet de choisir, parmi toute la population, les individus qui seront à
l’origine de la génération suivante. Il existe plusieurs méthodes. Les plus utilisées sont la
sélection par tournoi et la sélection proportionnelle.
• La sélection par tournoi :
25
Chapitre 1
j =1
- Le croisement
Cette étape est facultative pour les SE, et n’est pas toujours présente dans la génération de
nouveaux individus. La représentation par des nombres réels des individus rend les règles de
recombinaison simples.
• La recombinaison discrète
On choisit par sélection deux parents dans la population. Chaque composante du nouvel
individu est alors créée en choisissant aléatoirement la valeur de l’un des deux parents.
• La recombinaison intermédiaire
On choisit par sélection deux parents dans la population x1 et x 2 . Chaque composante du
nouvel individu x enfant est alors créée en faisant une combinaison linéaire des valeurs des deux
parents.
( )
x enfant i = xi1 + χ xi2 − xi1 (1.10)
Où χ est une variable aléatoire uniforme comprise entre 0 et 1, et fixée à 0.5 dans
beaucoup de cas.
26
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
- La mutation gaussienne
L’opérateur de mutation est toujours présent dans les SE. Il est le principal opérateur
d’exploration et est le garant de la recherche globale.
Le principe de la mutation gaussienne est d’ajouter un bruit gaussien centré sur les
variables que l’on souhaite faire muter.
xi = xi + N (0, σ ) (1.11)
Toute la difficulté réside dans le bon choix de la valeur de σ , σ étant la variance la loi
normale centrée N. La mutation auto-adaptative propose de faire évoluer cette grandeur avec
l’individu. Ainsi, les individus ayant une bonne valeur de σ effectueront des mutations
réussies et seront donc plus compétitifs. Cette grandeur σ évolue en suivant les mêmes
opérateurs que les variables d'optimisation. La mutation auto-adaptative a été proposée par
Rechenberg [Rec73] et Schwefel [Sch81]. Cette méthode présente l'avantage de ne pas ajouter
de paramètre de réglage de l'algorithme, et permet en plus d'assurer à la fois l'aspect global de
la recherche et une bonne efficacité de la recherche locale lors de la convergence vers la
solution.
27
Chapitre 1
coût seule ou valeur de la fonction coût conjuguée à son gradient. Dans ce dernier cas, la
continuité de la fonction coût est une nécessité.
Le choix de l'une ou l'autre des méthodes doit être adapté au problème à résoudre. Dans le
cadre de l'optimisation des assemblages par modélisation numérique, deux points primordiaux
doivent être retenus pour choisir la méthode d'optimisation :
• le temps de calcul de la fonction coût est long, car il fait appel à une résolution par
éléments finis non linéaire,
• le problème n'a pas de raison de présenter un seul minimum.
Les méthodes qui seront présentées par la suite tiendront compte de ces deux spécificités.
Dans un premier temps l'approche par plan d'expériences est présentée. Cette approche est
retenue car elle répond particulièrement bien à la problématique du coût de calcul : le nombre
d'évaluations de la fonction coût est parfaitement maitrisé.
Dans un second temps on se focalisera sur des méthodes d'optimisation dédiées à la
problématique du temps de calcul : stratégies d’évolution avec Métamodèle.
Historiquement la notion de PE date des années 1925, et a été initiée par Fisher et Yate
lors de leurs travaux en agronomie [Fis25]. La méthode gagne ensuite le monde industriel,
notamment grâce à Taguchi qui internationalise cette approche dans la fin des années 70
[Tag78] : il approfondit la méthode pour l'associer à la notion de qualité. Depuis, les PE ont
pris un essor considérable avec le développement de l'informatique et la puissance de calcul
qui l'accompagne.
28
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
Méthode d'Hotelling : dans cette méthode on mesure d'abord la somme des deux masses,
notée p1, puis la différence, notée p2.
ma + mb = p1 et ma − mb = p2
1 1 (1.12)
⇒ ma = [ p1 + p2 ] et mb = [ p1 − p2 ]
2 2
Si on applique le théorème de la variance à (1.12):
1 1 1
V( ma ) = V( p1 ) + V( p2 ) = σ 2 + σ 2 = σ 2 (1.13)
4 4 2
Où la variance est le carré de l'écart type de la mesure ( V( p ) = σ 2 ), on fait ici l'hypothèse
que l'erreur sur la mesure est de type gaussienne.
σ
L'erreur ainsi commise est de au lieu de σ , soit une amélioration de 30% de la
2
précision. On voit ici l'intérêt de faire varier tous les paramètres (les masses) à chaque
expérience (les pesées) : avec le même nombre d'expériences la précision est améliorée.
Comme le montre l'exemple précédent, une organisation rationnelle des expériences peut
réellement apporter une amélioration de la précision. Avant de détailler les différentes
méthodes de construction des plans d'expériences, présentons la démarche globale des PE.
Les 8 expériences et leurs réponses sont présentées dans le Tableau 1-2. Chacune de ces
expériences est réalisée avec une erreur σ = 2% au sens de la variance.
29
Chapitre 1
Ce type de plan est associé à une surface de réponse polynomiale Q1, dans notre cas :
y = m + E1x1 + E 2 x 2 + E 3 x 3 + I12 x1x 2 + I13 x1x 3 + I 23 x 2 x 3 + I123 x1x 2 x 3 (1.14)
Les effets sont directement à relier à la notion de sensibilité par rapport à un paramètre.
Les interactions caractérisent le couplage entre certains paramètres. Elles prennent tout leur
sens lorsque l'on étudie des réactions chimiques, où les interactions sont à l'origine de la
réaction.
Une représentation matricielle du plan est alors naturellement mise en place : Y le vecteur
réponse, E le vecteur des effets, et X la matrice d'expériences.
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y1 = 38 m
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y2 = 37 E1
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y3 = 26 E2
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y4 = 24 E3
X= ,Y= ,E= (1.15)
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y5 = 30 I12
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y6 = 28 I13
1 x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y = 19 I
7 23
1
x1 x2 x3 x1 x2 x1 x3 x2 x3 x1 x2 x3 y8 = 16 I123
Pour chaque ligne de la matrice X (1.15), x1, x2 et x3 prennent les valeurs +1 ou -1 suivant
l'essai réalisé (Tableau 1-2).
Le système s'écrit alors : Y = XE
Dans le cas d'un plan factoriel, la matrice d'expériences X a des propriétés spéciales
1
(matrices d'Hadamard), notamment X −1 = X t , où n est la taille de la matrice. On obtient
n
alors rapidement la valeur du vecteur d'effet :
30
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
m = 27.5
E1 = −1
E2 = −6
1 t E3 = −4
E= X Y⇒E= (1.16)
n I12 = −0.25
I13 = −0.25
I = 0.25
23
I =0
123
Cette dernière remarque fait apparaître que si l'on fait l'hypothèse a priori que cette
interaction est négligeable, un certain nombre d'essais peut être économisé. Ces types de plans
sont appelés plans factoriels fractionnaires.
Cette méthode, bien que coûteuse en nombre d'essais, est utilisée lorsque l'on traite peu de
paramètres. On peut citer le travail de Hana [Han06] sur le procédé de frettage : il réalise un
plan factoriel 23 pour modéliser le ratio d'énergie de glissement avec comme facteur la
pression de contact, la fréquence, et le déplacement. Dans ce cas les PE ont permis de
déterminer les conditions de frettage avec un minimum d'essais et de comparer deux états de
surface différents.
Où li est l'effet du facteur xi. On applique le même type de traitement que pour le plan
complet, les effets li sont alors calculés et comparés avec les effets Ei issus du plan complet
(Tableau 1-4).
31
Chapitre 1
Les résultats sont comparables, mais ils ont été obtenus avec moitié moins d'essais que
pour un plan complet. Cette réduction du nombre d'essais a forcément une contrepartie,
certaines approximations sont faites.
Détaillons le calcul de l3 :
1
l3 = [ − y2 − y3 + y5 + y8 ] (1.18)
4
Si l'on revient sur le plan complet :
1
E3 = [ − y1 − y2 − y3 − y4 + y5 + y6 + y7 + y8 ]
8
(1.19)
1
I12 = [ + y1 − y2 − y3 + y4 + y5 − y6 − y7 + y8 ]
8
On remarque alors : l3 = E3 + I12
l3 est en réalité la somme de l'effet E3 et de l'interaction I12, ces deux grandeurs sont dites
"aliasées". l3 est donc représentatif de l'effet du facteur x3 sous l'hypothèse que l'interaction I12
soit négligeable, ou inversement. La mise en place d'un plan fractionnaire doit donc être
précédée d'une phase de réflexion poussée sur les acteurs à prendre en compte, ainsi que sur
leurs interactions afin de décider lesquelles sont négligeables ou non.
Dans notre exemple, seuls les effets sont calculés, l'ensemble des interactions étant
considéré comme négligeable. Ce plan est dit de résolution I. Il est possible de réaliser des
plans de résolution supérieure. Le plan de résolution IV est très utilisé [Pil05, Chapitre 5] : les
facteurs sont "aliasés" avec les interactions entre 3 paramètres et les interactions d'ordre 2 sont
"aliasées" entre elles.
32
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
expériences déjà réalisées. L'alternative est donc de construire des plans optimaux, plans qui
donneront la meilleure précision sur les paramètres du modèle.
On se posera alors deux questions :
• que doit-on optimiser pour aboutir à un plan optimal ?
• suivant quel critère doit être réalisée cette optimisation ?
a) Précision du plan
Pour répondre à la première question, prenons le modèle linéaire suivant :
y = XE + e (1.20)
Ê est un estimateur de E. Cet estimateur par moindre carré a de très bonnes propriétés. On
démontre que sa matrice de variance-covariance est :
ˆ ( X t X )-1 σ 2
Cov(E)= (1.22)
où σ 2 = Var ( ε ) est la variance de l'approximation.
On arrive ainsi à exprimer la variance et la covariance des paramètres du modèle sans faire
intervenir les valeurs des réponses. On peut donc construire le plan d'expériences à priori. Un
PE optimal minimisera donc la valeur des termes de cette matrice.
b) Critères d'optimalité
Pour répondre à la deuxième question, citons les critères les plus utilisés.
• D-Optimalité : ce critère est basé sur la minimisation du déterminant de la matrice
X t X , que l'on peut assimiler à la variance globale du système, si l'on développe la
matrice de variance-covariance de Ê (1.21) :
( Xt X ) = det 1Xt X . Com ( Xt X )
−1 t
(1.23)
( )
où Com ( X ) est la co-matrice de X .
( ( ))
Le critère s'exprime alors π D = min det X t X , où U est l'ensemble des expériences
U
réalisables. Pour inclure la notion du nombre N d'expériences réalisées, on travaille avec le
déterminant normé.
π D = min
(
det X t X ) (1.24)
U N
Il est alors possible de déterminer le plan D-Optimal parmi un ensemble U de
N expériences possibles, pour un modèle donné. Ce type de critère est utilisé comme point de
départ pour la mise en place de plans à nombre d'expériences réduit [Kov07].
33
Chapitre 1
( (
p A = min trace X t X
U
)
-1
) (1.25)
34
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
Modélisation simplifiée
du procédé
(b)
Figure 1-7 : (a) Problème d'optimisation directe, (b) problème d'optimisation à deux niveaux
Il existe cependant des méthodes intrusives. Par exemple les méthodes de dérivation de
code : lors d'un calcul EF la valeur de la fonction coût est évaluée ainsi que son gradient. Ces
méthodes nécessitent un lourd travail dans le code source du logiciel EF. Pour plus de détails
sur ces méthodes on pourra se référer aux travaux de Forestier [For04].
Données du calcul EF
Paramètres d'entrée :
Variable d'optimisation Calcul EF Sortie :
Fonction coût
Boite noire
35
Chapitre 1
2.3.2 Métamodèle
De nombreuses techniques sont présentes dans la littérature pour réaliser cette
modélisation simplifiée. Le principe est de créer un opérateur mathématique, un métamodèle,
qui à partir d'une base de données de modélisation complète, propose une évaluation
approchée d'une solution non évaluée. Cet opérateur peut se représenter sous la forme d'un
schéma (Figure 1-9) ou d'une fonction (1.26).
fˆ ( x ) = Meta ( DBp, DBy, x ) (1.26)
Où fˆ est la valeur approchée de la fonction coût en x , DBp la base de données des jeux
de paramètres déjà évalués, DBy la base de données des fonctions coût associées.
Il existe une grande variété de métamodèles. Le Tableau 1-5, proposé par Bonte [Bon07],
permet d'avoir un aperçu de méthodes existantes ainsi que de leurs points forts et points
faibles. On pourra se référer à l'article de Gary Wang [Gar07] pour plus de détails sur les
différentes méthodes de construction de métamodèles.
36
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
37
Chapitre 1
Initialisation
x2
Calcul EF
D.E.
x1
Métamodèle Calcul EF
Méthode de minimisation
Non Oui
Convergence Solution
38
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
2.4 Bilan
Les méthodes présentées ici permettent toutes de résoudre un problème d'optimisation. Le
choix de l'une ou l'autre de ces méthodes se fait donc en analysant le type de problème que
l'on souhaite résoudre.
Les méthodes à direction de descente sont écartées, car elles sont mal adaptées aux
problèmes multi-extrema. A priori, rien ne nous permet de certifier que les problèmes
d'optimisation des assemblages et les problèmes d'identification de paramètres matériaux,
notamment de paramètres d'endommagement, sont mono-extrema. De plus ces méthodes à
direction de descente nécessitent que la fonction objectif soit dérivable ; la présence de
remaillage lors des simulations par éléments finis ne permet pas d'assurer une parfaite
continuité de la fonction coût.
Notre choix s'oriente donc vers des méthodes d'optimisation dites globales et ne
s'appuyant pas sur l'évaluation du gradient de la fonction objectif. Les méthodes issues des
algorithmes génétiques et par stratégie d'évolution semblent être de bonnes candidates pour
répondre à cette problématique. Elles présentent cependant l'inconvénient de nécessiter un très
grand nombre d'évaluations de la fonction coût. Ce nombre d'évaluations de la fonction coût
devient vite prohibitif lorsque le calcul de cette fonction coût fait appel à un calcul par
éléments finis.
Pour répondre à cette problématique supplémentaire de fonction objectif coûteuse en
temps de calcul, deux axes sont étudiés :
• La réduction du problème d'optimisation par une phase de sélection des
paramètres influents en s'appuyant par exemple sur les plans d'expériences.
• L'utilisation de méthodes d'optimisation assistées par métamodèle.
Finalement notre choix s'oriente sur l'exploitation de méthodes d'optimisation basées sur
les métamodèles par krigeage. En effet les algorithmes travaillant avec ce type de métamodèle
répondent aux différentes problématiques (recherche du minimum global, dérivabilité de la
fonction objectif non assurée, temps de calcul de la fonction coût important, méthode non
intrusive).
39
Chapitre 1
La base de données DBp est une liste de m jeux de paramètres. Chaque jeu de
paramètres étant un vecteur de dimension n ( n représente la dimension du problème
d'optimisation). La base de données des réponses DBy est, quant à elle, un vecteur de
m valeurs.
Cette base de données doit être initialisée, on fait à ce stade l'hypothèse qu'elle contient
déjà m points.
Pour construire l'estimateur ŷ du processus gaussien Fx , m et σ 2 sont évalués par la
méthode des moindres carrés :
1 t .C −1 . y
m = t −1 (1.30)
1 .C . 1
σ=
(y − 1 .m ) .C .(y − 1 .m )
t −1
(1.31)
n
Cette matrice de corrélation C est construite à l'aide des fonctions de corrélation. L'une
des fonctions la plus utilisée pour l'analyse d'expériences numériques est une fonction de
forme exponentielle. Une formulation en est donnée à l'équation (1.33) :
(
corr (x, x′,θ ,τ ) = exp − l ( x, x′)
τ
) (1.33)
40
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
avec
2
n
x′ − x
l ( x, x′ ) = ∑ i i (1.34)
i =1 θi
La fonction de corrélation exponentielle est retenue car elle permet d'obtenir une surface
de réponse continue et dérivable si τ =2, et continue mais non dérivable si 1 ≤ τ < 2 . Les
procédés d'assemblage s'appuient sur des phénomènes physiques continus, il est donc légitime
de choisir une forme de surface de réponse continue.
(a) Fonction gaussienne (ττ=2), équation (b) Fonction exponentielle non dérivable
(1.33) (ττ=1), équation (1.33)
41
Chapitre 1
3.1.2 Calibration
Pour illustrer les aspects liés à la calibration, un exemple est utilisé. Une surface de
réponse 1D par krigeage est construite à partir de 4 points maîtres.
Si θ est trop faible la surface obtenue n'a que peu d'intérêt, la seule information disponible
est alors la valeur de la tendance m (équation(1.30)). Chaque point maître est lié à ce plan
42
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
médian par un pic très local (Figure 1-13a, pointillés noirs). Ce plan médian de valeur
m apparaît nettement sur la Figure 1-13a.
Si la valeur de θ est bien calibrée la surface de réponse est beaucoup plus riche en
informations (Figure 1-13a, trait continu bleu).
Le paramètre τ permet de choisir entre une surface de réponse dérivable (Figure 1-13b,
trait bleu continu) et une surface non dérivable (Figure 1-13b, pointillés noirs). Cette propriété
de dérivabilité est directement liée à la forme des fonctions de corrélation. La fonction de
corrélation est dérivable en zéro uniquement pour τ =2. La surface de réponse associée est
donc dérivable, aux points maîtres, uniquement pour τ =2.
b) Méthodes de calibration
Il existe principalement deux méthodes pour réaliser cette calibration. La méthode de
validations croisées (Cross Validation, CV) et la méthode de l'estimateur du maximum de
vraisemblance (Maximum Likelihood Estimation, MLE). Ces méthodes s'appuient sur la base
de données des points maîtres. Elles ont pour objectif de minimiser un critère en jouant sur les
paramètres de calibration θ et τ .
43
Chapitre 1
Cette méthode a pour principal avantage de ne faire appel qu'au déterminant de la matrice
de corrélation C . La construction complète de l'opérateur de krigeage n'est pas nécessaire. La
difficulté résidant dans la forme de la fonction à maximiser. Cette fonction est multimodale,
des méthodes adaptées doivent donc être utilisées pour la maximiser. Ulmer et al. [Ulm03]
proposent d'utiliser une méthode évolutionnaire, Laurenceau et al. [Lau08] proposent quant à
eux d'utiliser une méthode à direction de descente avec un point de départ spécifique, évalué à
partir de l'analyse du phénomène modélisé.
L'objectif de toutes ces méthodes est de proposer un point maître supplémentaire à ajouter
dans la base de données. L'ajout successif de points maîtres aboutit à la résolution du
problème de minimisation.
La recherche de ce point a pour effet de préciser un minimum déjà détecté. Cette phase est
appelée exploitation. Sur la
l'approximation par krigeage est représentée en trait continu bleu, le point pm +1 proposé
par cette méthode est représenté sur la
par le point pm +1 vert, le minimum détecté dans cette zone sera donc connu de manière
plus précise.
44
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
−1
pm+1 : Exploration
y
−2
−3 pm+1 : Exploitation
−4
Points maitres
−5 Approximation y^
Variance σ^
−6
−1 −0.5 0 0.5 1
x
L'exploitation et l'exploration sont deux aspects qui doivent être pris en compte pour
construire un algorithme d'optimisation global. La combinaison de (1.39) et (1.40) permet de
créer, par addition, un critère multi-objectifs. Ce critère est nommé Lower confidence Bound
(LB) [Emm06].
pm +1 = xmin | xmin = min ( yˆ ( x ) − ωσˆ ( x ) ) , ω ∈ [1,3] (1.41)
x∈ X
Ce critère introduit un coefficient de pondération ω qui doit être fixé par l'utilisateur. Plus
ω est grand plus l'effet d'exploration est accentué. La Figure 1-15 (trait fin noir) représente le
critère LB ( ω =2). Dans cet exemple le point pm +1 minimisant le critère se trouve en x = 0.3 .
Dans la zone se situant entre −1 ≤ x ≤ −0.5 , l'approximation est considérée comme fiable, la
part exploitation du critère est donc faible, l'exploration est donc favorisée.
Le critère LB ajoute un coefficient supplémentaire à déterminer ω , la multiplication de ces
coefficients de "réglage" étant toujours dommageable pour la simplicité d'utilisation des
méthodes.
45
Chapitre 1
Il existe d'autres méthodes qui ne font pas appel directement à une somme pondérée de
ŷ ( x ) et σˆ ( x ) . Le critère d'amélioration espérée (Expected Improvement, EI) proposé par
Jones [Jon98] en est un exemple.
y − yˆ ( x ) y − yˆ ( x )
EI ( x ) = ( ymin − yˆ ( x ) ) Φ min + σˆ ( x ) ϕ min
σˆ ( x )
(1.42)
σˆ ( x )
Où ymin est le minimum des réponses de la base de données ( ymin = min ( yi ) ), Φ la
i∈[1, m ]
Sur l'exemple Figure 1-15 (trait rouge pointillé) l'EI présente deux pics. L'un en x = 0.3
correspond à l'aspect exploration, et l'autre en x = −0.25 correspond à l'aspect exploitation.
Dans ce cas les aspects globaux et locaux de la recherche de l'optimum sont bien présents. Ce
critère a l'intérêt d'exploiter les propriétés de l'opérateur de krigeage sans ajouter de paramètre
de "réglage" supplémentaire.
Exploration
Exploitation Exploration
Figure 1-15 : Critère Lower confidence Bounds (LB) et Expected Improvement (EI)
46
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
L'obtention des coordonnées d'un point maître supplémentaire se fait toujours via la
minimisation d'un critère (ou maximisation suivant les cas). Cette minimisation n'est jamais
triviale. L'exemple présenté en fil rouge de ce paragraphe est très simple (1 variable
d'optimisation, 4 points maîtres) et pourtant la forme des critères à minimiser est complexe.
Les fonctions à minimiser sont multi-extrema. L'EI Figure 1-15 en est un exemple, l'un des
extrema est même relativement étroit. Ces minimisations font appel à des méthodes
d'optimisation bien choisies. Les méthodes par stratégie d'évolution (paragraphe 2.1) sont
particulièrement bien adaptées. La résolution de ce problème d'optimisation annexe n'est pas
bloquant en terme de temps. En effet, cette résolution est bien moins coûteuse en temps de
calcul que l'évaluation de la fonction coût liée à un point maître obtenue par un calcul
éléments finis.
Dans les travaux de Jones, la maximisation de l'Expected Improvement est réalisée par
l'algorithme branch-and-bound algorithm [Jon98], qui est spécialement adapté à la forme de
la fonction EI ( x ) .
Le schéma de la Figure 1-16 propose une vue du fonctionnement de l'algorithme
d'optimisation globale. Dans ce schéma la maximisation de EI ( x ) est réalisée par un
algorithme à stratégie d'évolution.
47
Chapitre 1
Initialisation
Plan
d’expériences
P (1…m)
P (1…m)
Calcul E.F.
Exploration de la
surface de réponse
f(1…m) Optimisation par S.E.
Opérateurs
génétiques
EI Nouvelle
⇒ EI ( x) population
oui non
P_min m+1
fin
48
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
• l'utilisation d'un critère capable de proposer Npar points maîtres par itération,
• la création d’un système pseudo-itératif où, à chaque pseudo-itération, un point
maître est proposé, ce point maître étant associé à une valeur de fonction coût
fictive et temporaire.
La première alternative est encore peu présente dans la littérature. Ginsbourger et al.
[Gin10], proposent un critère sur la base de l'Expected Improvement, nommé q-EI, où q est le
nombre de points maîtres extraits à l'aide de ce critère. La formulation de ce critère reste assez
complexe et n'est calculée analytiquement que pour q=2. Pour des valeurs supérieures de q, le
critère est évalué par méthode de Monte-Carlo.
Une deuxième alternative est décrite dans le Tableau 1-7. Cette approche consiste à
construire tout d'abord un métamodèle s'appuyant sur la base de données.
Ce métamodèle est ensuite exploité pour extraire un nouveau point maître (Tableau 1-7,
ligne 9). A ce stade le nouveau point maître n'est pas évalué directement via l'appel au calcul
en "boite noire". Ce point maître devient un point maître fictif : la valeur de sa fonction coût
associée est fixée de manière arbitraire (Tableau 1-7, ligne 11).
Un nouveau métamodèle est alors construit sur la base de données des points maîtres
augmentée du point maître fictif (Tableau 1-7, ligne 8). Un nouveau point maître fictif est
alors créé. Cette opération est répétée autant de fois que l'on souhaite de points maîtres fictifs.
Une fois ce nombre atteint, les fonctions coût associées à l'ensemble des points maîtres
fictifs sont calculées de manière exacte en faisant appel au calcul EF en "boite noire" (Tableau
1-7, ligne 14). Ces Npar calculs peuvent alors être exécutés simultanément sur différents
processeurs.
Après avoir évalué ces Npar points maîtres fictifs, ils deviennent des points maîtres inclus
dans la base de données (Tableau 1-7, ligne 15). En itérant ainsi, les points maîtres évalués
sont choisis par "paquets" de Npar nouveaux points par pseudo-itération.
49
Chapitre 1
50
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
f * ( x ) = yˆ ( x ) (1.44)
L'enrichissement par KB est illustré Figure 1-17. L'exemple traite de l'ajout de deux
points. L'étape 1, Figure 1-17a, est de déterminer le nouveau point par maximisation de l'EI
(équation(1.42)). La valeur de la fonction coût associée à ce point maître est fixée à la valeur
prédite par l'opérateur de krigeage (Figure 1-17b, point rouge).
Un nouveau métamodèle est alors construit sur la base des points maîtres et du point
maître fictif. L'ajout du point maître fictif ne modifie pas la valeur de l'estimateur (Figure
1-17b, trait bleu), en revanche l'estimation de l'erreur est modifiée. Sa valeur est nulle pour
tous les points maîtres, donc automatiquement au point maître fictif. Le critère de l'EI est
donc modifié. Le nouveau point déterminé par la maximisation de l'EI est donc dans une autre
zone (Figure 1-17b, point vert).
La dernière étape est de calculer la valeur des fonctions coût des nouveaux points maîtres
fictifs ainsi rajoutés (Figure 1-17c, point rouge).
51
Chapitre 1
(c) Etape 3
52
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
f * ( x ) = yL (1.45)
Plusieurs valeurs de yL peuvent être choisies : le minimum ( yL = min ( yi ) ), le maximum
i∈[1, m ]
( yL = max ( yi ) ) ou la moyenne ( yL = moy ( yi ) ) des valeurs des fonctions coût des points
i∈[1, m ] i∈[1, m ]
maîtres.
Plus yL est grand et plus la recherche sera globale (pour un problème de minimisation).
Cette méthode a pour principal risque d'introduire temporairement un point maître fictif
aberrant dans la base de données. Ce problème est illustré Figure 1-18. A la première étape
(Figure 1-18a, point vert) le point maximisant l'EI est sélectionné. A l'étape 2 (Figure 1-18b,
point rouge) une valeur fictive de fonction coût est fixée suivant la méthode CL-max
(équation(1.45)). Or cette valeur est très différente de la valeur du point maître le plus proche.
Le métamodèle étant interpolant, cette variation très brusque se répercute par de fortes
oscillations du métamodèle (Figure 1-18b, trait bleu), le métamodèle devient inexploitable. Il
est cependant à noter que dès le calcul des fonctions réelles effectué, le métamodèle redevient
valide et exploitable (Figure 1-18c).
53
Chapitre 1
(c) Etape 3
54
Optimisation dans le cadre de la modélisation par éléments finis des procédés d'assemblage
5 Conclusion
Dans ce chapitre, après une présentation succincte de la problématique industrielle liée
aux procédés d’assemblage par déformation plastique, les principes généraux sur
l'optimisation ont été présentés, en se focalisant sur l'aspect global ou local de la recherche.
Dans un deuxième temps, la problématique du temps d'évaluation de la fonction coût a été
abordée. En effet, dans le cadre de l'optimisation des assemblages, le test d'une solution
d'assemblage fait appel à la résolution d'un problème éléments finis non linéaire. Le temps
d’obtention de la valeur de la fonction coût associée est donc grand.
L'approche par plans d'expériences est une première réponse à cette problématique. Sans
toutefois permettre de construire un processus d'optimisation, elle permet de maîtriser le
nombre d’évaluations de la fonction coût. Elle peut également s’avérer utile pour évaluer la
sensibilité de la fonction coût aux paramètres d’optimisation, ou encore pour construire une
première base de données de points maîtres.
Une seconde réponse est l'exploitation des méthodes basées sur les métamodèles. Ces
méthodes introduisent un niveau de modélisation simplifié qui est utilisé pour résoudre un
problème d'optimisation approché plus rapide. Le métamodèle par krigeage a été plus
particulièrement présenté, ainsi que l'algorithme EGO.
La dernière réponse proposée est d'utiliser la parallélisation pour accélérer le temps de
résolution d’un problème d’optimisation. Dans le cadre de l’algorithme EGO avec utilisation
du métamodèle par krigeage, l'objectif est de pouvoir évaluer la fonction coût de plusieurs
points maîtres de manière simultanée.
Le chapitre suivant présente la plateforme MOOPI (MOdular software dedicated to
Optimization and Parameters Identification) développée dans le cadre de cette thèse. Nous
nous intéresserons plus particulièrement à l’implémentation des algorithmes présentés dans ce
chapitre.
55
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Chapitre 2
MOOPI un outil d’optimisation dédié aux calculs éléments
finis
MOOPI
Jeu de paramètres
à évaluer
Fonction coût
Logiciel de calcul
"Boite Noire"
57
Chapitre 2
58
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Remarque : ces deux étapes peuvent être remplacées par un calcul analytique
"rapide" de la fonction coût associée aux paramètres d'entrée. Cette substitution est
faite lors des tests de développement afin de réduire le temps imputé au calcul de la
fonction coût.
59
Chapitre 2
Entrées Sorties
OBSERVABLES
PARAMETRES
Optimisation
Paramètres
Couche 2 optimaux
60
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Ce schéma (Figure 2-2) fait apparaître en premier lieu une base commune à tous les
objectifs : la couche 0, l'exécution du calcul et l'évaluation de la fonction coût. En effet, tous
les types d'analyses (analyse de sensibilité, optimisation, ou analyse inverse) peuvent être vus
comme des méthodes de gestion d'un ensemble de calculs à réaliser.
Cette dernière observation pousse à faire le choix de placer cet ensemble de calculs à
réaliser au centre de l'organigramme de la plateforme.
61
Chapitre 2
1.4.3 L'explorateur
L'explorateur a comme rôle d'analyser la base de données des calculs. Grâce à cette
analyse l'explorateur propose de nouveaux jeux de paramètres à tester pour atteindre l'objectif.
L’explorateur pourra par exemple être une méthode d’optimisation implémentée dans la
plateforme MOOPI.
Ces trois blocs permettent de définir les trois classes de base de l'architecture orientée
objet de la plateforme MOOPI. Une quatrième classe permet de définir l'environnement
global de la plateforme et d'organiser les actions des trois blocs.
L'objectif est d'exploiter la puissance de calcul disponible pour limiter, au plus, le temps
pour atteindre l'objectif fixé. Ces gestionnaires de calculs prennent tout leur sens lorsqu'ils
62
MOOPI un outil d’optimisation dédié aux calculs éléments finis
sont associés à des méthodes d'optimisation parallèles telles que les méthodes présentées dans
le Chapitre 1.
Calculs à lancer
63
Chapitre 2
64
MOOPI un outil d’optimisation dédié aux calculs éléments finis
L'utilisation des classes virtuelles est, en plus d'un guide pour le programmeur, un moyen
de sécuriser le code, puisque l'ajout d'une nouvelle méthode ne justifie en rien la modification
des autres classes.
Figure 2-6 : Communication par mots clés - extrait de fichier de mise en données Forge
L’optimisation peut ainsi être réalisée sur tout type de variable utilisée dans le fichier de
données : paramètres matériaux, conditions aux limites, conditions d’interfaces, etc.
L’optimisation de forme, quant à elle, ne peut fonctionner que si la géométrie optimisée est
paramétrée.
Un travail un peu plus conséquent est nécessaire pour construire la fonction coût. Dans le
schéma de communication de la plateforme, la valeur de la fonction coût doit être écrite sous
forme d'un scalaire dans un fichier texte. Il peut donc être nécessaire dans certains cas
d'ajouter un module supplémentaire de post-traitement, qui traite les résultats du logiciel de
calcul, puis construit la fonction coût.
65
Chapitre 2
66
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Base de données
des calculs
Explorateur EGO
Gestion du
critère d'arrêt
Algorithme d'optimisation par
Stratégie d'Evolution
Rôle :
• rechercher les coordonnées du
maximum de l'EI
Base de données
des calculs
67
Chapitre 2
Points maîtres
Krigeage
Construction Construction
Calibration
yˆ ( x )
Evaluation en un point x Evaluation σˆ ( x )
EI (x )
68
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Krigeage
Stratégie d'évolution
Population initiale
Construction
Opérateur génétique
(sélection, croisement, mutation)
Convergence
Coordonnées
du minimum
69
Chapitre 2
Nombre de nouveaux
Base de données des points maîtres
calculs Explorateur EGO
Points maîtres
Krigeage
Points maîtres
virtuels
S.E.
Parallélisation
Affectation d'une valeur virtuelle
de la fonction coût au nouveau
point maître
Non
Base de données
des calculs
Nb de nvx
Nouveaux points points maîtres 1 nouveau point
maîtres atteints maître
Oui
70
MOOPI un outil d’optimisation dédié aux calculs éléments finis
71
Chapitre 2
non
NbCalcul oui
Fin
<NbCalculMax
Nouveau point
maître non
oui
EI>0. 1% Itération suivante
non
Ce critère, où l'on autorise la recherche de points par exploration pure, permet de limiter
les risques de convergence prématurée vers un minimum local, convergence prématurée qui
peut avoir lieu si la base de données initiale est trop pauvre.
72
MOOPI un outil d’optimisation dédié aux calculs éléments finis
73
Chapitre 2
Minimum 15
50
10
10
0
0
global
10 0
10
10
0
10 50
50
10
y
50
100
Minimum
50
local 5 10
10
50
10
10
10
0
−5 0 0 5 10
x
74
MOOPI un outil d’optimisation dédié aux calculs éléments finis
2 0.5
0.5
0.01
Minimum 1.5
1
0. 0
1 01
1
0.0 0.0
y
0.5
1
0.01
0.5
0.5
0
1
0 0.5 1 1.5 2
x
Figure 2-13 : Iso valeurs de la fonction de Rosenbrock 2D
3 2
H 3,4 = −∑ α i exp −∑ Aij ( x j − Pij ) + 3.86278,
4
i =1 j =1
1 3 10 30 6890 1170 2673 (2.3)
1.2 0.1 10 35 4699 4387 7470
α = , A = , P = 10−4
3 3 10 30 1091 8732 5547
3.2 0.1 10 35 381 5743 8828
Le fonction de Hartmann 3D est définie pour 0 < xi =1,2,3 < 1 . Cette fonction présente 4
minima locaux et un minimum global situé en x = (0.11, 0.55, 0.85) et valant 0. Une
illustration de cette fonction est proposée Figure 2-14 (seulement 2 minima apparaissent sur
cette vue).
75
Chapitre 2
i =1 j =1
1 10 3 17 3.05 1.7 8
1.2 0.05 10 17 0.1 8 14
α = , A = , (2.4)
3 3 3.5 1.7 10 17 8
3.2 17 8 0.5 10 0.1 14
1312 1696 5569 124 8283 5886
4135 8307 3736 1004 9991
−4 2329
P = 10
2348 1451 3522 2883 3047 6650
4047 8828 8732 5743 1091 381
La fonction de Hartmann 6D est définie pour 0 < xi =1..6 < 1 . Cette fonction présente 6
minima locaux et un minimum global situé en x = (0.20, 0.15, 0.47, 0.27, 0.31, 0.65) et valant 0.
Les deux fonctions de Hartmann permettent de tester les deux aspects de la recherche pour
des dimensions plus élevées : l'exploration (la fonction est multi extrema), et l'exploitation (le
gradient proche de l'optimum est faible, ce mauvais conditionnement apparait sur la Figure
2-14).
76
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Point de départ
1
6 2
5 3
Figure 2-15 : Test de l’algorithme BFGS sur la fonction de Branin-Hoo (carré vert : minimum
global, ronds rouges : itérations de la méthode BFGS)
15
10
50
10
Minimum
0
0
0
global
0
10
10
0 10
10 50
50 Minima
10
50
locaux
y
50
100
10
5
10
50
10
10
10
0
0
−5 0 5 10
x
Figure 2-16 : Test de l'algorithme EGO sur la fonction de Branin-Hoo (carré vert : minimum
global, étoiles rouges : points maîtres initiaux, ronds rouges : itérations de la méthode EGO)
Ces deux premiers tests montrent que la plateforme MOOPI répond à l'objectif général de
résolution d'un problème de minimisation présentant plusieurs minima.
77
Chapitre 2
78
MOOPI un outil d’optimisation dédié aux calculs éléments finis
79
Chapitre 2
Déroulement séquentiel
Solution
Solution
Solution
Itérations / temps
0 1 2 3 4 5
Un pas d'optimisation
Remarque : La méthode de calcul du speed-up fait l'hypothèse que tous les calculs de la
fonction coût ont la même durée (Figure 2-19, déroulement parallèle parfait). Or lorsque le
calcul de la fonction coût fait appel à un calcul éléments finis, les temps d'exécution peuvent
être variables le speed-up réel sera alors moins bon (Figure 2-19, déroulement parallèle réel).
Ces tests de speed-up sont réalisés sur les fonctions coût tests de référence. Le speed-up
est calculé pour p appartenant à la plage {1, 2, 4,8,16} .
3.4.2 Résultats
Les résultats présentés montrent une étude comparative entre deux types d'enrichissement
virtuel de la base de données. Le premier nommé Kriging Believer (KB) et le second Constant
Liar avec comme valeur constante le minimum des valeurs des points maîtres (CL-min). Ces
deux types d'enrichissement virtuels sont définis au Chapitre 1, dans la Partie 4.3.
Les résultats des fonctions testées sont donnés en Figure 2-20. En pointillé rouge est
représenté le speed-up théorique S p = p , en trait plein croix noire le speed-up avec
enrichissement CL-min et en trait plein cercle noir le speed-up avec enrichissement KB.
80
MOOPI un outil d’optimisation dédié aux calculs éléments finis
(a) (b)
(c) (d)
(e)
81
Chapitre 2
82
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Population initiale
Évaluation de la population
Krigeage : Calcul EF :
Evaluation Evaluation
approché exacte
Opérateur génétique
(sélection, croisement, mutation)
Krigeage : Calcul EF :
Evaluation Evaluation
approchée exacte
Convergence
Coordonnées
du minimum
83
Chapitre 2
(a) base de données initiale de 2 points - (b) base de données initiale de 2 points (ronds
convergence noirs) - répartition des points testés
(c) base de données initiale de 4 points - (d) base de données initiale de 4 points (ronds
convergence noirs) - répartition des points testés
(e) base de données initiale de 9 points - (f) base de données initiale de 9 points (ronds
convergence noirs) - répartition des points testés
Figure 2-22 : Comparaison EGO-p / SE-Meta sur la fonction de Branin-Hoo - 1 calcul à chaque
itération
84
MOOPI un outil d’optimisation dédié aux calculs éléments finis
Les tests menés dans les mêmes conditions, en considérant 4 calculs exactes de la fonction
coût par itération ou par génération pour SE-Meta (Figure 2-23) mènent à la même conclusion
de forte dépendance de l'algorithme SE-Meta à la base de données initiale. Il est tout de même
à noter que l'utilisation de 4 calculs de la fonction par itération, pour SE-Meta, permet une
meilleure exploration de l'espace de recherche (Figure 2-23 b et d) sans pour autant assurer la
détection du minimum global.
(a) base de données initiale de 4 points - (b) base de données initiale de 4 points -
convergence répartition des points testés
(c) base de données initiale de 9 points - (d) base de données initiale de 9 points -
convergence répartition des points testés
Figure 2-23 : Comparaison EGO-p / SE-Meta sur la fonction de Branin-Hoo - 4 calculs à chaque
itération
85
Chapitre 2
(a) base de données initiale de 9 points - (b) base de données initiale de 9 points -
convergence répartition des points testés
Figure 2-24 : Comparaison EGO-p / SE-Meta sur la fonction de Branin-Hoo - 8 calculs à chaque
itération
Pour assurer un bon comportement de SE-Meta (sur ce problème), il faut donc évaluer 8
individus de manière exacte à chaque génération. La fonction coût de ces 8 points peut alors
être évaluée de manière simultanée.
Pour les tests suivants, on se place donc dans cette configuration, qui est la seule assurant
la résolution du problème d'optimisation : à chaque génération de SE-Meta, 8 calculs exacts
de la fonction doivent être effectués.
Pour évaluer les performances on définit le temps apparent pour résoudre le problème de
minimisation, où p est le nombre de processeurs utilisés :
nbCalcul
t* = (eq. 2.6)
p
Si on dispose uniquement d'une machine monoprocesseur, l'algorithme EGO-p est deux
fois plus rapide que SE-Meta (Figure 2-25a). En effet pour assurer le bon comportement de
SE-Meta, les 40 calculs (sélectionnés 8 par 8) nécessaires sont lancés de manière totalement
séquentielle. Alors que EGO-p (avec p=1) ne nécessite que 19 calculs pour converger (calcul
sélection 1 par 1).
Si l'on dispose de la puissance pour évaluer 4 calculs simultanément alors les
performances de EGO-p (avec p=4) sont 1.5 fois supérieures à SE-Meta (Figure 2-25b).
Si, en revanche, on dispose de la puissance de calcul pour évaluer ces 8 fonctions coût
simultanément alors les performances des deux algorithmes sont comparables en termes de
vitesse de convergence et de globalité de la recherche (Figure 2-25c).
86
MOOPI un outil d’optimisation dédié aux calculs éléments finis
87
Chapitre 2
Là encore la résolution est plus rapide avec l'algorithme EGO-p. On peut tout de même
noter la convergence de SE-Meta vers le minimum global.
88
MOOPI un outil d’optimisation dédié aux calculs éléments finis
4 Conclusion
Dans ce chapitre ont été présentés, dans un premier temps, les développements réalisés
pour la mise en place de la plateforme MOOPI. Le développement de la plateforme, et
l'écriture des plus de 5000 lignes de code en C++, ont permis de répondre aux deux objectifs
fixés au début de ce chapitre : la mise en place d'un outil d'optimisation dédié à l'exploitation
de calculs éléments finis, ainsi que la mise en place d'un environnement capable d'intégrer des
méthodes d'optimisation diverses et de gérer plusieurs logiciels de calcul scientifique.
Contrairement à l’algorithme BFGS, l’algorithme EGO-p permet d’identifier le minimum
global, ce qui nous sera très utile dans le cadre de l’identification de paramètres.
Dans un deuxième temps, l'algorithme d'optimisation EGO-p, présenté au Chapitre 1, a été
testé. Le résultat majeur obtenu sur les fonctions tests est une accélération significative pour
un nombre de processus parallèles de l'ordre du nombre de dimensions du problème
d'optimisation. Cette accélération permet d'aboutir à la solution d'un problème d'optimisation
plus rapidement. L'objectif principal étant de proposer un outil capable de réaliser une
optimisation sur la base de calculs par éléments finis lourds en terme de temps calcul. De
plus, l'un des points forts de la méthode EGO-p est le faible nombre de paramètres à calibrer
pour effectuer l'optimisation.
La plateforme développée au cours de cette thèse est vouée à être enrichie. Les axes
principaux d'enrichissement sont relatifs aux méthodes d'optimisation proposées. La structure
a été pensée pour accueillir des méthodes d'optimisation plus avancées telles que les méthodes
multi-objectifs, ainsi que des méthodes d'optimisation sous contraintes. L'algorithme EGO-p
peut lui aussi faire l'objet d'extensions, en intégrant par exemple le traitement d'un bruit
éventuel sur la fonction coût. Cette extension peut être envisagée en enrichissant le
métamodèle par krigeage de l'effet dit de "pépite", qui permet d'associer à chaque point maître
une incertitude sur sa valeur : le métamodèle n'est alors plus interpolant ([Kle09] [Bee03]).
Les développements présentés dans ce chapitre permettent de répondre à l'un des objectifs
majeur de cette thèse. Objectifs qui étaient de développer un outil capable de répondre aux
problématiques d'optimisation du projet MONA LISA mais aussi plus largement aux
problématiques d'optimisation au sein du laboratoire ainsi que de capitaliser les différents
développements réalisés.
Les chapitres suivants s'appuieront sur la plateforme MOOPI pour réaliser de
l'optimisation d'assemblages mécaniques ainsi que l'identification de paramètres matériaux
par analyse inverse. D'autres exemples d'application sont rassemblés en Annexe de ce
manuscrit.
89
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Chapitre 3
Optimisation numérique d’un point d’assemblage :
application au procédé de clinchage
Simulation Simulation
Pose de l’assemblage Histoire Tenue mécanique en service
Mécanique
Optimisation
91
Chapitre 3
Poinçon
Tôle supérieure
Tôle Inférieure
Bouterolle
Les points forts de ce type d'assemblage sont multiples : la mise en œuvre est simple et
économique, il n'y a pas de problèmes de compatibilité entre les matériaux comme c'est le cas
pour le soudage par point. Ce type d'assemblage connait une forte croissance dans le domaine
de l'automobile. Mais son utilisation reste cependant limitée par manque de procédures
standards et d'outils de dimensionnement [Jom07].
La modélisation numérique de ce procédé reste peu développée, comme indiqué dans la
revue de litterature de He [He10] sur la modélisation éléments finis du clinchage. On peut tout
de même citer les travaux de Hamel [Ham00] sur la modélisation de la pose et les travaux de
Jomaâ [Jom07] sur la modélisation de la tenue mécanique tenant compte du passé
thermomécanique issu de la mise en place de l'assemblage. Il existe cependant une
bibliographie plus large sur la modélisation d'autres technologies d'assemblage. On peut par
exemple citer les articles de Bouchard et al. [Bou08] et de Porcaro et al. [Pro06] sur le
procédé de rivetage auto-poinçonneur ou encore l'article de Chen et al. [Che11] sur le procédé
de rivetage.
Dans ce chapitre on présente tout d'abord la modélisation d'un assemblage par clinchage,
incluant une validation de la modélisation de la phase de pose ainsi que de la tenue mécanique
par rapport à des essais expérimentaux. Sur la base de cette modélisation validée, une analyse
de sensibilité de la tenue mécanique par rapport aux paramètres du procédé est réalisée. Cette
analyse de sensibilité permet ensuite de sélectionner les paramètres les plus influents et de
réaliser l'optimisation de la tenue de l'assemblage vis-à-vis de ces paramètres influents,
optimisation réalisée grâce à la plateforme MOOPI.
92
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
• Le comportement élastique :
σ = λ tr ( ε e ) I d + 2 µε e (3.2)
où λ et µ sont les coefficients de Lamé, et I d la matrice identité. Les relations entre les
coefficients de Lamé et le module d'Young E et le coefficient de Poisson ν sont rappelées :
E Eν
µ= ,λ = (3.3)
2 (1 + ν ) (1 +ν )(1 − 2ν )
93
Chapitre 3
Germination
Ces trois phases sont identifiables sur une courbe de traction, elles sont représentées en
Figure 3-4. Cette courbe illustre l'impact de l'apparition de ces cavités sur le comportement du
matériau. Dès la phase de croissance (segment CD, Figure 3-4) la courbe présente un
adoucissement dû à la croissance des micro-cavités. Cet adoucissement s'amplifie lors de la
phase de coalescence (segment DE, Figure 3-4), la courbe force/déplacement s'écarte
nettement de la courbe hypothétique ne prenant pas en compte cette évolution (segment CG
pointillé, Figure 3-4).
Remarque : L'adoucissement progressif observé sur une courbe force/déplacement relative
à un essai de traction (Figure 3-4) est dû à l'action conjointe de la croissance de
l'endommagement et de la striction de l'éprouvette.
94
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Dans le cadre du projet MONA LISA, il a été choisi de travailler avec un modèle
phénoménologique pour représenter ce phénomène. Ce type de modèle semble bien adapté au
contexte industriel du projet et a déjà montré des résultats très satisfaisants dans des travaux
précédents [Bou08, Bou10].
Dans ce type d'approche phénoménologique, le phénomène est modélisé par une variable
d'endommagement macroscopique w qui peut être définie comme le rapport suivant :
S − Sɶ S D
wn = = (3.6)
S S
95
Chapitre 3
appliquée, ramenée à la section qui résiste réellement aux efforts appliqués. Elle est définie de
la manière suivante :
σ
σɶ = (3.7)
(1 − w )
La contrainte effective fait apparaitre la contrainte σ , qui est la contrainte dans le matériau
dépourvu d'endommagement. Cette contrainte est issue, dans notre cas, de la résolution d'un
problème mécanique élasto-plastique décrit précédemment.
La formulation de la loi d'évolution de l'endommagement présentée ici introduit plusieurs
apports, tels que la limite de triaxialité en compression et la décomposition du tenseur des
contraintes en parties positive et négative.
L'objet ici n'est pas de proposer une description détaillée du modèle d'endommagement
ductile, mais simplement de donner les clés nécessaires pour la mise en place de la
modélisation des assemblages mécaniques dans un premier temps puis dans un deuxième
temps pour l'identification des paramètres du modèle. Pour plus de détails, le lecteur pourra se
référer aux travaux de thèse de R. El Khaoulani [Kha10] menés en parallèle dans le cadre du
projet MONA LISA, ainsi qu’aux travaux de L. Bourgeon [Bou10].
λ pl Y b
− si 0 ≤ Tx et si ε > ε d
1 − w S 0
λ pl Y
b
1
wɺ = − si − ≤ Tx < 0 et si ε > ε d (3.8)
1 − hw S0 3
1
0 si Tx ≤ − ou si ε ≤ ε d
3
wɺ = 0 et w = 1 si w ≥ wc
96
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
−1 (1 +ν ) σ : σ −ν tr (σ ) 2
Y=
2 E (1 − w )
2
+ +
(3.10)
h
2 (
− 1 + ν ) σ : σ −ν tr (σ ) 2
2 E (1 − w )
− −
L'allure de cette loi d'évolution est représentée Figure 3-6. Sur cette figure, apparaissent
deux des paramètres du modèle : wc le seuil d'endommagement critique, et ε d la déformation
seuil de déclenchement. Les paramètres caractérisant l'endommagement du matériau sont
donc b, S0, ε d et wc.
wc
0 ε
εd
97
Chapitre 3
98
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Cette variable non locale est définie sur tout le domaine d'étude Ω et fait intervenir une
fonction de pondération ϕ . La fonction de pondération choisie est de forme gaussienne :
1 x− y 2
ϕ ( x, y ) = exp − (3.13)
( 2π ) lc 3 2lc
32 2
99
Chapitre 3
w ( x ) − lc ∆w ( x ) = w ( x )
2
(3.14)
∇w ( x ) .n = 0 sur ∂Ω
Figure 3-8 : Iso valeurs d'endommagement – approche non locale par gradient implicite
[Kha10]
100
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Dans cette partie, on définit tout d'abord le modèle éléments finis utilisé pour la simulation
des deux phases. On s'attardera sur l'importance de la conservation de l'histoire
thermomécanique du matériau pour obtenir une bonne corrélation entre les résultats issus de
la simulation et ceux issus d'essais expérimentaux.
Une fois cette modélisation mise en place et validée, les aspects liés à l'optimisation du
procédé seront abordés : le choix des paramètres à optimiser, les grandeurs observables, ainsi
que la mise en place d'une chaine de simulations automatisée.
Paramètres Valeur
σy 30 MPa
K 464MPa
n 0.32
b 1
S0 2.25 MPa
εd 0.17
wc 0.8
h 0.1
Tableau 3-1 : Paramètres matériau pour la modélisation du clinchage
a) Paramètres du procédé
En plus des deux tôles, cette phase met en jeu 3 outils ( Figure 3-9) :
• Le poinçon : cet outil est piloté en déplacement avec une vitesse de 1 mm/s, sa
course est stoppée sur un critère de force.
• La matrice (ou bouterolle): cet outil est fixe.
101
Chapitre 3
• Le serre flan (ou dévêtisseur) : cet outil est un outil flottant, il est piloté en force, la
force Fz est évaluée suivant la formule :
Fz = F0 + K d ∆U (3.15)
où F0 est une force constante, K d la raideur du ressort compris dans le dévêtisseur,
et ∆U le déplacement relatif entre le poinçon et le serre flan [Rey07].
Les caractéristiques de cette loi sont évaluées sur la base des essais expérimentaux, la
raideur est évaluée à 690kN/m, la force initiale à 2000N.
Z
Fz
Tôle 1
Tôle 2
Matrice
102
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
(a) (b)
Pour modéliser l'éventuelle rupture des tôles une méthode par kill-element est utilisée.
Lorsqu'un élément atteint l'endommagement critique, cet élément est supprimé (ce
phénomène n'apparait pas sur la Figure 3-10).
c) Validation expérimentale
La configuration présentée ici a fait l'objet d'une campagne d'essais en collaboration avec
le CETIM dans le cadre du projet MONA LISA. La mise en place de l'assemblage a été
réalisée sur des installations industrielles spécialement instrumentées pour mesurer la force et
le déplacement du vérin pilotant le poinçon au cours du procédé de pose.
Les courbes présentées ici sont les courbes force/déplacement ramenées au poinçon, c'est-
à-dire que l'on a retiré le déplacement dû à la déformation élastique et aux jeux du dispositif,
et la contribution en force absorbée par le dévêtisseur. 36 poses de la même configuration ont
été réalisées. La Figure 3-11 présente les deux courbes extrêmes enregistrées.
Une seconde grandeur est mesurée à la fin du procédé de clinchage : la côte X. La côte X
est l'épaisseur finale cumulée des deux tôles en tête du poinçon après la pose de l'assemblage
103
Chapitre 3
(Figure 3-14). Les mesures réalisées sur les 36 éprouvettes sont illustrées Figure 3-12. La
distribution de cette grandeur est en accord avec les tolérances industrielles qui sont de 15%.
1.4
Figure 3-11 : Effort de pose expérimental et numérique d'un assemblage par clinchage
Figure 3-12 : Distribution des mesures de côte X sur les essais de pose de point clinché
104
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
dans la première phase ( 0 < d < 1.4mm ) dû à une rigidité trop forte des tôles en flexion,
cependant cet écart reste minime. Cet écart peut éventuellement être expliqué par une
surestimation du frottement ou encore par une surestimation de l'écrouissage à faible valeur
de déformation.
Outre la bonne corrélation de la courbe d'effort de pose, la côte X obtenue par simulation
est de 0.47 mm. Cette grandeur est en accord avec la plage des mesures expérimentales
présentées Figure 3-12.
Les résultats présentés ici retranscrivent bien les points caractéristiques de la courbe de
pose :
• L'inflexion à d = 1.4mm correspond au contact de la tôle inférieure avec le fond de
la matrice (Figure 3-13a).
• La phase de compression (Figure 3-13b). Lors de cette phase l'accroissement de
l'endommagement est le plus fort (la tôle supérieure est en traction au dessus du
rayon de congé du poinçon).
• La phase de remplissage de la gorge de la matrice (Figure 3-13c) où la force
augmente rapidement (d>2.8mm). Lors de cette dernière phase, l'endommagement
augmente principalement sur la tôle inférieure au fond de la gorge de la matrice, et
le S se forme entre les deux tôles ; S qui permet de créer l'assemblage des deux
tôles.
105
Chapitre 3
Collet
106
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Outil 2
Z
Cote X Tôle 2
Tôle 1
Outil 1
La modélisation ainsi obtenue est présentée Figure 3-15. Le calcul est réalisé sur un PC
(Intel(R) Core(TM)2 CPU 6700 @ 2.66GHz, 2Go Ram) en un temps de 24 min.
(a)
(b)
107
Chapitre 3
(c)
Figure 3-15 : Modélisation axisymétrique de la tenue mécanique dans la direction Z d'un point
clinché (a : état initial, b : après rupture au collet, c : champs d'endommagement juste avant
rupture)
Le résultat obtenu est présenté Figure 3-15, et la courbe de l'effort d'arrachement associée
est tracée Figure 3-16. Ce graphe montre que le calcul de la force est perturbé (Figure 3-16,
points bleus). Ces perturbations sont dues à la gestion du contact dans la zone de glissement
entre les deux tôles. Pour exploiter cette courbe, et en extraire une valeur représentative du
maximum, il est nécessaire de lisser ces perturbations.
Pour effectuer ce lissage on utilise un filtre de convolution gaussien, décrit équation (3.17)
.
i N
F
lissée = ∑
j =− N
i− j
Fbrute .g j
1 1 − x2
2 (3.17)
g =j
e , x = [ − N , − N + 1,..., N − 1, N ]
∑i
g i 2π
i i
Où Fbrute et Flissée sont respectivement les mesures de force brutes et lissées, gi les
coefficients de pondération et 2N+1 le nombre de points de mesure utilisé pour le filtrage.
Ce filtre réalise une moyenne pondérée glissante des données brutes. N est fixé à 10 pour
le filtrage de la force, c'est-à-dire que chaque point est réévalué en réalisant la moyenne
pondérée des 20 mesures entourant ce point. La pondération de chaque point est donnée par la
loi gaussienne centrée réduite g j (équation(3.17)). La courbe filtrée est présentée Figure 3-16
(trait noir). La force maximale ainsi évaluée s'affranchit des valeurs perturbées par la gestion
du contact. La mise en place de ce filtre est primordiale pour évaluer la force d'arrachement et
utiliser cette grandeur comme observable à maximiser dans le processus d'optimisation.
108
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
109
Chapitre 3
Réinitialisation de l'endommagement
Historique complet
Sans historique
Les courbes de l'effort d'arrachement sont présentées Figure 3-17 (courbes lissées).
L'impact du transfert de l'historique thermomécanique est très net : si l'on ne transfère pas cet
historique, la résistance mécanique est sous évaluée d'un facteur 2. L'exportation, ou non, de
la variable d'endommagement a un impact plus faible sur la force d'arrachement. Elle est
surestimée de 20% si l'on néglige l'historique relatif à l'endommagement. De plus, la prise en
compte de l'endommagement dans le passé thermomécanique permet de bien prédire le type
de rupture du point d'assemblage. La Figure 3-18a présente le mode de rupture, en tenant
compte de tout le passé thermomécanique. La rupture du point a lieu au collet, c'est le mode
de rupture qui a été observé expérimentalement. En revanche les Figure 3-18b et c, qui ne
tiennent pas compte de la totalité du passé thermomécanique, présentent un désassemblage et
non une rupture.
110
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
(a) transfert de tous les champs (b) transfert de tous les champs sauf
l'endommagement
Pour valider cette modélisation, des essais de tenue mécanique ont été réalisés. Ils ont été
réalisés sur une machine de traction équipée d'un montage de type ARCAN. Ces essais nous
permettent d'obtenir la force d'arrachement du point d'assemblage dans la direction du
poinçon (axe Z). L'essai a été répété 4 fois, et les forces d'arrachement sont présentées dans le
Tableau 3-2. Ces résultats présentent une dispersion relativement importante. Cela est
notamment dû à la variabilité des propriétés matériau des tôles, ainsi qu'à la variabilité déjà
introduite par le procédé de pose. La force d'arrachement de 737N prédite par la simulation
est en bon accord avec les résultats expérimentaux.
111
Chapitre 3
Révolution
112
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Le maillage utilisé dans la zone proche du point d'assemblable est le maillage issu de la
pose, il est simplement obtenu par révolution. Ce maillage est présenté Figure 3-20.
Six essais de tenue mécanique dans cette configuration ont été réalisés à l'aide d'un
montage ARCAN. Les valeurs d'arrachement obtenues sont présentées dans le Tableau 3-3.
Là encore les résultats présentent une forte dispersion.
113
Chapitre 3
Les paramètres à optimiser sont relatifs aux outils : dimensions et pilotage. Ils sont listés
dans le Tableau 3-4 et illustrés Figure 3-21.
Paramètres Notation
Poinçon
Force maximum de pose Fmax
Rayon Rp
Angle sous tête W1
Rayon de congé Rcp
Angle de dépouille W2
Matrice
Profondeur Pm
Largeur gorge Lm
Dévêtisseur
Raideur Kd
Tableau 3-4 : Paramètres à optimiser pour le clinchage
114
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
Z w2
Lm
Z
Rcp Pm
w1
Rp
(b) Matrice
(a) Poinçon
Les paramètres variables sont donc au nombre de huit. Le nombre de paramètres liés à la
matrice est limité car la génération automatique du profil de l'outil s'avère relativement
complexe sans outil de CAO dédié.
MOOPI
Jeu de
paramètres à
évaluer
Fonction coût
Logiciel de calcul
"Boite Noire"
Cette chaîne de simulations est donc mise en place, elle est décrite sur l’organigramme
présenté en Figure 3-23. Ce chaînage fait se succéder plusieurs opérations réalisées avec
115
Chapitre 3
116
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
8 Paramètres d'entrée
Pose
Préparation du calcul
Incorporation des outils
modifiés
(GLpre® en mode console)
Exécution du calcul
(Forge2D-2009®)
Tenue
mécanique Préparation du calcul
Calage des outils
(GLpre® en mode console)
Exécution du calcul
(Forge2D-2009®)
Post-process
Récupération de la force
d'arrachement (filtrage)
(Matlab®)
117
Chapitre 3
3 Analyse de sensibilité
La première partie de ce chapitre a permis de mettre en place une modélisation fiable du
procédé de clinchage, et de la tenue mécanique. Huit paramètres procédé ont été considérés
comme variables potentielles. Il n'est cependant pas envisageable, et même non pertinent, de
réaliser l'optimisation de la tenue mécanique de ces 8 paramètres. Une pré-étude est donc
réalisée pour discriminer les paramètres les plus influents. Cette phase est appelée phase de
screening [Bon07].
On fait le choix ici de réaliser une analyse de sensibilité en évaluant le gradient de la force
d'arrachement par rapport aux 8 paramètres procédé (Tableau 3-4). Les valeurs prises par
chacun des paramètres sont présentées dans le Tableau 3-5.
Les 9 calculs (pose et tenue mécanique décrits Figure 3-23) sont lancés à l'aide de la
plateforme MOOPI. Les forces d'arrachement correspondant à chaque calcul sont présentées
Tableau 3-5, ainsi que la valeur du gradient.
Ces résultats sont présentés sous forme d'un diagramme de Pareto (Figure 3-24). Ce
diagramme présente la sensibilité relative à chacun des paramètres, dans l'ordre décroissant de
la valeur du gradient (histogramme). Il présente aussi la courbe du cumul de chacun de ces
gradients, exprimée en pourcentage (trait noir).
118
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
0.35 100
90
0.3
80
0.25 70
60
% cumulé
0.2 df/dp
Df/Dp
%cumulé 50
0.15
40
0.1 30
20
0.05
10
0 0
Pm Rp Lm Rc W2 Fmax W1 Kd
p
Cette analyse permet de déterminer les paramètres les plus influents et donc ceux qu'il est
le plus pertinent d'optimiser. Deux paramètres géométriques sont donc retenus pour la suite de
l’étude d’optimisation : la profondeur de la matrice (Pm) et le diamètre du poinçon (Rp).
119
Chapitre 3
Valeur
Plage de
Valeur nominale Optimale
Variation
identifiée
Rp (mm) 1.9 [1.6 2.2] 1.96
W1(°) 5 -
Rc (mm) 0.3 -
W2(°) 2 -
Pm (mm) 1.2 [-0.3 0.6] 0.16
Lm (mm) 1.1 -
Fmax (kN) 2.03 -
Kd(N/m) 690000 -
Tableau 3-6 : Plage de variation des paramètres pour l'optimisation et résultats
4.2 Résultats
Les résultats de l'optimisation sont présentés Tableau 3-6, Figure 3-25 et Figure 3-26.
Tout d'abord, nous pouvons noter que la méthode converge rapidement vers une zone où
la force d'arrachement est supérieure à 800N. Cette valeur est obtenue en 4 itérations, soit 8
calculs du modèle direct, en plus des 3 calculs de la phase d'initialisation (Figure 3-25). Les
itérations suivantes présentent ensuite une amélioration de la fonction coût beaucoup plus
lente.
Après 14 itérations, le critère de l'Expected Improvement indique la convergence de la
méthode : la surface de réponse est connue de manière fiable. La solution optimale est
obtenue à la 8ieme itération. Les paramètres optimaux sont présentés Tableau 3-6, et la force
d'arrachement associée est de 840N. L'optimisation des deux paramètres Rp et Pm permet
donc un gain de 13.5% de la tenue mécanique.
120
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
900
800
700
600
500
F(N)
400
300
200
Fc -−Initialisation
Fc initilsiation
100 Fc − iteration
Fc − min
0
0 5 10 15
Iterations
121
Chapitre 3
Surface de Reponse
200
400
400
0.5 600
0
60
0.4
800
0.3
0
80
6 00
Pm
0.2 800
600
400
0.1
0
40
0 200
600
400
0
−0.1 200
Maximum detecté
0
−0.2 Config. ref.
Iso−F(N)
L'observation des coupes des points après assemblage (Figure 3-27) permet d'expliquer le
gain obtenu en terme de tenue mécanique. La solution optimisée propose une formation du S
moins prononcée, 0.11 mm contre 0.17 mm pour la configuration de référence. Cette tendance
va pourtant à l'encontre d'une augmentation de la tenue mécanique. En revanche la tenue
mécanique est directement impactée par le taux d'endommagement induit par la phase de pose
dans la zone du collet : la configuration optimisée réduit de 60% à 30% l'endommagement
dans cette zone. Cette réduction d'un facteur 2 explique le gain de tenue mécanique, malgré
l'atténuation du S.
Le rayon du poinçon est légèrement supérieur dans la configuration optimisée par rapport
à la configuration de référence, et pourtant la formation du S est moins prononcée. Cela
provient de la profondeur plus importante de la matrice qui crée un espace d’écoulement de la
matière au cours du clinchage, et diminue ainsi l’endommagement dans la zone du collet.
0.17 0.11
0.18 0.22
122
Optimisation numérique d’un point d’assemblage : application au procédé de clinchage
5 Conclusion
Dans ce chapitre une procédure d'optimisation de la tenue à l’arrachement d'un point
d'assemblage par clinchage a été présentée. Cette optimisation a été réalisée en faisant varier
certains paramètres géométriques des outils de clinchage. Les points clés de cette procédure
sont :
• La validation de la modélisation du procédé et de la tenue mécanique du point
clinché.
• La mise en place d'un chaînage automatique de ces deux modélisations, en
conservant l'historique thermomécanique du matériau.
• L'étude de sensibilité dans l'objectif de discriminer les paramètres procédé les plus
influents.
• L'optimisation de la tenue mécanique en elle-même.
L'optimisation a mené à un gain de 13.5% sur la tenue mécanique dans la direction Z, gain
obtenu en modifiant les deux paramètres les plus influents : le rayon du poinçon et la
profondeur de la matrice. Cette optimisation réalisée sur la tenue mécanique dans l'axe Z a
permis, de plus, d'augmenter la tenue mécanique du point en cisaillement de 42.1%. Une
123
Chapitre 3
124
Identification de paramètres matériaux par analyse inverse
Chapitre 4
Identification de paramètres matériaux par analyse
inverse
Les chapitres précédents ont permis de mettre en place un outil d'optimisation dédié au
calcul par éléments finis. Cet outil a pour objectif principal d'optimiser la tenue mécanique
des assemblages en tenant compte des phases de pose et de tenue en service, tout en
conservant l'histoire thermomécanique des différents matériaux entre ces deux phases. Pour
réaliser cette optimisation (Figure 4-1), il est important d'avoir une bonne connaissance du
comportement mécanique des matériaux mis en jeu.
Dans ce chapitre, la plateforme développée est utilisée pour l'identification par analyse
inverse des paramètres des lois de comportement retenues. La majeure partie du chapitre sera
consacrée à l'identification des paramètres de la loi d'endommagement de Lemaitre [Lem01]
couplée à la loi de comportement élasto-plastique du matériau. Contrairement à
l’identification de paramètres élastiques ou élasto-plastiques en petites déformations,
l’identification de paramètres d’endommagement en grandes déformations plastiques ne
permet pas l’utilisation de méthodes analytiques à partir d’essais normalisés. Il existe des
méthodes expérimentales, directes (coupes et observations du taux de porosité, tomographie
…) ou indirectes (mesures de module d’Young, de densité, de résistivité électrique …),
permettant d’évaluer l’état d’endommagement dans un matériau au cours de la déformation.
Ces mesures restent cependant très délicates à mettre en place, et la corrélation entre les
observations et les paramètres d’un modèle (phénoménologique dans notre cas) sont
discutables. C’est pourquoi on ne s’intéressera ici qu’à l’approche basée sur l’analyse inverse.
125
Chapitre 4
Loi de comportement
Simulation Simulation
Pose de l’assemblage Histoire Tenue mécanique en service
Mécanique
Optimisation
Observable
Le schéma Figure 4-2 peut aussi se formuler sous la forme d'une fonction :
Obs = F ( P ) (4.1)
126
Identification de paramètres matériaux par analyse inverse
Le modèle direct, noté F, est donc capable d'associer à des paramètres d'entrée P une
observable de sortie Obs.
Dans certain cas, les paramètres d'entrée sont inconnus et seule l'observable de sortie est
connue. On souhaite alors connaitre les paramètres d'entrée qui ont permis la génération de
cette observable. Formellement on souhaite inverser la fonction F (définie en (4.1) ) :
P = F −1 ( Obs ) (4.2)
Cette inversion est impossible à expliciter dans un grand nombre de cas. On cherche alors
à trouver le jeu de paramètres P qui satisfait l'équation (4.3)
F ( P ) − Obs = 0 (4.3)
fc ( P ) = F ( P ) − Obs (4.4)
Remarques :
127
Chapitre 4
MOOPI
Jeu de paramètres
à évaluer
Fonction coût
Logiciel de calcul
"Boîte Noire"
Cette "boîte noire" est décrite en Figure 4-4. La fonction coût est construite par
comparaison de deux observables, l'une numérique issue du modèle direct, l'autre
expérimentale.
Boîte noire
Observable
expérimentale
Observable numérique
Paramètres Modèle Déformations
modèle Contraintes Fonction Coût
Direct Courbe force/déplacement
…
128
Identification de paramètres matériaux par analyse inverse
La nature des observables, ainsi que la construction de la fonction coût seront détaillées
ultérieurement.
2 Le modèle direct
Il est aussi important de bien différencier une grandeur observable, d'une grandeur pilote
de la simulation. Par exemple lors d'un essai de traction piloté en déplacement, le déplacement
est la grandeur pilote et la force une grandeur observable. Inversement lors d'un essai de
129
Chapitre 4
La spécification de chacun des paramètres du Tableau 4-1, ainsi que des paramètres
numériques, permet de définir le modèle direct.
Avant d'exploiter des informations très riches issues de méthodes de mesure avancées, on
souhaite exploiter les informations disponibles avec un essai de traction classique. Ce type
d'essai est relativement simple à mettre en place, et est très répandu dans le domaine
industriel. Cette première approche permettra de plus de bien expliciter le module d’analyse
inverse, et de montrer comment la fonction coût peut être exploitée pour mettre en évidence la
présence de minima locaux.
130
Identification de paramètres matériaux par analyse inverse
( )
n
fc ( P ) = ∑ Obsinum ( P ) − Obsiexp
2
(4.6)
i =1
où n est le nombre de points de mesure, Obsinum et Obsiexp sont respectivement les valeurs
des mesures des observables numériques et expérimentales, et P le jeu de paramètres testé.
Cette formulation est très largement utilisée dans la résolution de problèmes inverses ; elle
est par exemple utilisée par Springmann et Kuna [Spr03] pour identifier les paramètres du
modèle d'endommagement de Rousselier.
∑ ( Obs ( P ) − Obs )
n 2
num exp
i i
fc ( P ) = i =1
(4.7)
∑ ( Obs )
n
exp 2
i
i =1
Cette formulation adimensionnelle est présente dans les travaux de Harth et al. [Har04].
Cette formulation permet notamment à Hart de construire une fonction sur la base de
plusieurs observables issues de différents essais en réalisant directement une somme de ses
fonctions coût.
131
Chapitre 4
∑ ( Obs ( P ) − Obs )
n 2
num exp
i i
fc ( P ) =
i =1
(4.8)
n 2
( ) ( )
2 n
min ∑ Obsi exp
, ∑ Obsi ( P )
num
i =1 i =1
Cette formulation permet de ne plus avoir une limite finie de la fonction coût lorsque
l'observable numérique tend vers 0 (Figure 4-5, trait pointillé vert).
fc
Obsnum
Obsexp
Figure 4-5 : Fonction coût normalisée - illustration du cas où Obsnum tend vers 0, trait plein bleu
formulation (4.7), trait pointillé vert formulation (4.8)
∫ ( Obs ( P, x ) − Obs ( x ) )
2
num exp
dx
fc ( P ) = x
(4.9)
( ) (
min ∫ Obs exp ( x ) dx, ∫ Obs num ( P, x ) dx )
2 2
x x
132
Identification de paramètres matériaux par analyse inverse
( )
n
( ) ( )
n n
(4.10)
∑ ∑ ( )
2 2
min Obsi
exp
∆xi , Obsi
num
P ∆xi
i =1 i =1
∆xi = xi − xi −1
Cette formulation permet d'exploiter des observables issues d'une résolution éléments finis
à pas de temps variable, et des observables expérimentales issues de capteurs à taux
d’échantillonnage variable.
On pourra par la suite ajouter une fonction ω ( xi ) de pondération si l'on souhaite donner
un poids plus important à une partie de l'observable. La fonction coût utilisée dans la suite de
ce chapitre prend donc la forme suivante :
( )
n
ω ( xi ) . Obsinum ( P ) − Obsiexp ∆xi
2
∑
fc ( P ) = i =1
n
( ) ( )
n
(4.11)
∑ ( ) exp 2
i ∑ ω ( xi ) .Obsi ( P ) ∆xi
2
i =1 i =1
∆xi = xi − xi −1
La fonction coût présentée ci-dessus (4.11) fait l'hypothèse que les mesures des
observables sont coïncidentes, c'est-à-dire que pour une courbe force/déplacement chaque
mesure de force Fi, expérimentale et numérique, est réalisée pour un même déplacement di.
Or, dans la plupart des cas, cette hypothèse n'est pas vérifiée. Pour obtenir ces points
coïncidents, une interpolation linéaire de l'observable expérimentale est réalisée.
a) Problématique
La Figure 4-6 présente, schématiquement, les observables expérimentale et numérique
d'un essai de traction. Le calcul de la fonction coût dans la zone 1 (Figure 4-6) ne pose pas de
problème particulier. En revanche, le calcul de la fonction coût dans la zone 2 (Figure 4-6) est
plus délicat. Proche de la zone de rupture, pour un déplacement donné, l'une des deux
observables n'est plus disponible. C'est-à-dire que la rupture a eu lieu pour deux allongements
différents. Cette zone est pourtant une zone primordiale pour l'identification des paramètres
de la loi d'endommagement ; il est donc nécessaire de bien évaluer l'écart entre les deux
observables dans cette zone.
133
Chapitre 4
F Obsexp
Obsnum
Zone 1 Zone 2
Afin d'illustrer les différentes évolutions du calcul de la fonction coût, on considère les
paramètres élasto-plastiques connus, et on trace la surface de réponse de la fonction coût
suivant deux paramètres du modèle d’endommagement : S0 et ε d . La fonction coût est
évaluée sur la base d'une observable expérimentale virtuelle, c'est-à-dire issue directement du
modèle direct pour un jeu de paramètres donné (ce jeu de paramètres étant censé représenter
la solution optimale de l’analyse inverse). Ce type de test sera présenté plus en détail dans la
partie 3.3 de ce même chapitre. Les plages de recherche des paramètres sont normalisées pour
varier entre -1 et 1. Les paramètres optimaux sont situés au centre du domaine en (0,0), le
domaine est représenté dans l'échelle réduite [-1,1].
b) Première solution
La première solution envisagée est de calculer la fonction coût sur la fenêtre où les deux
observables sont disponibles, c'est-à-dire la zone 1 de la Figure 4-6. La surface de réponse
obtenue est présentée Figure 4-7.
134
Identification de paramètres matériaux par analyse inverse
Fonction coût V1
0.03
0.02
Fc
0.01
0
1
−1
0 0
−1 1 Epsd*
S0*
• Pour des valeurs de S0* faibles c'est-à-dire proche de -1, la fonction coût décroît
vers un autre minimum que celui situé au centre du domaine. Cette zone est
représentative d'un cas pathologique de cette méthode de calcul (Figure 4-8). Si la
rupture est très rapide, que la phase d'adoucissement est quasi inexistante, le calcul
de la fonction coût est réalisé uniquement sur la première partie de la courbe. Dans
cette partie les deux observables sont très proches (on considère ici que les
paramètres élasto-plastiques sont connus). Un écart apparaît uniquement lors de la
rupture, mais cet écart à un impact très faible dans le calcul de la fonction coût. En
effet, en plus d'être très bref, cet écart est pondéré par un pas d'échantillonnage
∆x (4.11) très faible.
135
Chapitre 4
F Obsexp
Obsnum
Rupture
Zone de calcul de Fc
Figure 4-8 : Description schématique des observables force/déplacement - cas d'une phase
d'adoucissement quasi inexistante
• Pour des valeurs de ε d * fortes, la fonction coût reste constante lorsque S0* varie.
Cette zone correspond à un autre cas pathologique (illustré Figure 4-9). Les deux
observables numériques, représentées Figure 4-9, ont un allongement à la rupture
différent, mais cette différence se situe hors de la zone de calcul de la fonction
coût, les valeurs des fonctions coût associées à ces deux observables numériques
sont donc identiques. Cela explique l'apparition d’une zone où la fonction coût est
constante sur la surface de réponse Figure 4-7.
F Obsexp
Obsnum1 Dernier point pour calcul
Obsnum2
Zone de calcul de Fc
136
Identification de paramètres matériaux par analyse inverse
• La surface de réponse est bruitée. Ce bruit est bien visible dans la zone où S0* est
grand. Ce bruit est principalement du à la position du dernier point de l'observable
numérique utilisé pour le calcul de la fonction coût (4.11). Ce problème est illustré
Figure 4-9. Si le point repéré par le cercle rouge est utilisé ou non pour le calcul de
la fonction coût (suivant son abscisse di) la valeur de celle-ci sera légèrement
perturbée, ce qui aboutit à une réponse bruitée.
Cette première approche n'est pas acceptable ; les trois problèmes soulevés ici rendent très
difficile la résolution du problème de minimisation.
F Obsexp
Obsnum
∆F
Points de mesure
ajoutés
d
dmax
Compléter les mesures permet de résoudre les 3 problèmes observés sur la Figure 4-7. En
effet, aucune des trois pathologies n'est visible sur la surface de réponse obtenue en
complétant les mesures, Figure 4-11.
137
Chapitre 4
Fonction coût V2
0.6
0.4
Fc
0.2
0
1
−1
0 0
−1 1 Epsd*
S0*
La fonction coût, calculée en complétant les mesures par des points de mesure de force
nulle, prend des valeurs très importantes si l'on s'éloigne de la solution optimale. En effet, on
ajoute des termes de l'ordre de ∆F 2 (Figure 4-10) dans la fonction coût. Une forte différence
de gradient apparaît sur la surface de réponse (Figure 4-11), forte différence qui induit un
problème mal conditionné. Ce qui pénalise le processus de minimisation.
Pour limiter ce phénomène on propose de compléter les mesures par des points, non plus
de force nulle, mais de la valeur du dernier point de mesure de l'observable expérimentale Ffin
(Figure 4-12). Cette dernière modification permet d'aboutir à une surface de réponse (Figure
4-13) avec des différences de gradient plus faibles, et donc plus simple à exploiter par la
procédure de minimisation.
F Obsexp
Obsnum
Ffin
Points de mesure
ajoutés
138
Identification de paramètres matériaux par analyse inverse
Fonction coût V3
0.6
0.4
Fc
0.2
0
1
−1
0 0
−1 1 Epsd*
S0*
3.3.1 Méthodologie
Afin de réaliser cette étude de sensibilité, on travaille sur la base d'un essai virtuel. C'est-à-
dire que l'observable numérique est générée par une exécution du modèle direct avec un jeu
de paramètres matériau donné, qui sera donc le jeu de paramètres à identifier. Travailler avec
une observable virtuelle permet d’une part de connaître la solution du problème
d'identification et donc d'analyser l'efficacité de la méthode. Cela permet d’autre part de
s'affranchir des incertitudes expérimentales et numériques, puisque les deux observables sont
générées dans les mêmes conditions, en terme de paramètres numériques. Cela permet enfin
d’utiliser le même maillage et le même pas de temps pour s'affranchir d'éventuelles
dépendances numériques et ainsi se concentrer sur la méthode d'identification par analyse
inverse.
139
Chapitre 4
Paramètres Jeu 1
E 69000 MPa
ν 0.3
σy 46 MPa
K 430MPa
n 0.34
b 1
S0 0.7 MPa
wc 0.8
εd 0.16
lc 0.06 mm
h 0.2
Tableau 4-2 : Jeu de paramètres de référence pour la génération d'un essai virtuel
205
105
76
25
12.5
R=20
(a)
Capteurs de déplacement
utilisés pour mesurer l'allongement
(b)
Figure 4-14 : Eprouvette de traction normalisée pour produit plat - (a) dimensions (mm) - (b)
maillage éléments finis
Pour ce test, un quart de l'éprouvette est modélisé, un maillage fin dans la zone de
localisation de l'endommagement est utilisé (Figure 4-14b). La courbe force/déplacement
ainsi obtenue est tracée Figure 4-15. Le déplacement est mesuré à l'aide de deux capteurs de
140
Identification de paramètres matériaux par analyse inverse
déplacement (Figure 4-14b), capteurs qui permettent de recréer une mesure d'allongement
comparable à une mesure expérimentale par un extensomètre. L'allongement est mesuré sur le
centre de la zone utile et sa longueur initiale est de 50mm.
3000
2500
F(N) 2000
1500
1000
500
0
0 2 4 6 8 10 12
d(mm)
Les surfaces de réponse ainsi obtenues sont présentées Figure 4-16 à Figure 4-21.
141
Chapitre 4
Déplacement / Force
2.2 1.8
0.18
0.07
2 0.16
1.6
0.14 0.06
1.8
0.12 0.05
1.6 1.4
0.1
b
b
0.04
1.4 0.08 1.2
0.03
1.2 0.06
0.04 1 0.02
1
0.02 0.01
0.8 0.8
1 2 3 0.4 0.6 0.8 1
S0 S0
Figure 4-16 : Surface de réponse - observable : Figure 4-17 : Surface de réponse - observable :
force – paramètres étudiés : b/S0 force – paramètres étudiés : b/S0- détail
Déplacement / Force Déplacement / Force
2.2 0.1 6
0.09 0.08
2
5
0.08 0.07
1.8
0.07 0.06
4
1.6 0.06 0.05
b
b
0.05
1.4 3 0.04
0.04
0.03
1.2 0.03
2 0.02
0.02
1
0.01 0.01
0.8 1
0 0.1 0.2 0 0.05 0.1 0.15 0.2
Epsd Epsd
Figure 4-18 : Surface de réponse - observable : Figure 4-19 : Surface de réponse - observable :
force – paramètres étudiés : b/ ε d force – paramètres étudiés : b/ ε d - plage élargie
du paramètre b
142
Identification de paramètres matériaux par analyse inverse
Déplacement / Force
3 1.2
0.16 0.12
2.5 0.14
1 0.1
0.12
2 0.08
0.1 0.8
S0
S0
1.5 0.08 0.06
0.6
1 0.06 0.04
0.04 0.4
0.5 0.02
0.02
0.2 0
0 0.1 0.2 0.12 0.14 0.16 0.18 0.2
Epsd Epsd
Figure 4-20 : Surface de réponse - observable : Figure 4-21 : Surface de réponse - observable :
force – paramètres étudiés : S0/
εd εd
force – paramètres étudiés : S0/ - détail
• Figure 4-16, Figure 4-17 : le paramètre S0 présente une forte sensibilité sur ces
différentes surfaces de réponse. La valeur optimum de S0 est facilement détectable
pour une valeur du paramètre b donnée. En revanche, le paramètre b a une
sensibilité très faible sur la fonction coût, cette sensibilité est cependant visible sur
la Figure 4-17 où le domaine d'étude est réduit autour du jeu de paramètres
optimum. Cette surface ne montre pas de couplage important entre les paramètres
b et S0 dans la zone de l'optimum.
• Figure 4-18, Figure 4-19 : cette surface de réponse présente une forte sensibilité
par rapport au paramètre ε d lorsque b est faible (b<1.5). Dans cette même zone, la
sensibilité par rapport au paramètre b est faible. La surface de réponse Figure 4-19,
où la plage de variation de b a été agrandie, fait apparaître un minimum local
(b=2.5, ε d =0.14). La sensibilité autour de ce minimum local est relativement
faible. Deux séries de calculs supplémentaires ont été réalisées pour analyser ce
minimum local. La Figure 4-22 présente l'évolution de la fonction coût en fonction
de b, pour deux valeurs de ε d . Ces deux cas sont matérialisés par les traits
pointillés sur la Figure 4-18 et la Figure 4-19. Les deux minima apparaissent. La
valeur de la fonction coût relative à ce minimum local est plus élevée que celle
relative au minimum global mais l'écart est faible.
143
Chapitre 4
0.04
Fc Fc
0.035
Epsd=0.16
0.03
Epsd=0.14
0.025
0.02
0.015
0.01
0.005
bb
0
0.5 1 1.5 2 2.5 3
Figure 4-22 : Evolution de la fonction coût en fonction de b pour deux valeurs de εd données
• Figure 4-20, Figure 4-21 : la surface de réponse présente une sensibilité forte par
rapport aux paramètres S0 et ε d . Sur les deux figures, des minima locaux
apparaissent mais ils sont le fruit du tracé des iso-valeurs de la fonction coût (voir
remarque suivante). La zone de l'optimum, la vallée, est de forme courbe, et ne
semble pas particulièrement orientée suivant un axe principal. L'identification
couplée de ces deux paramètres est donc obligatoire lorsque l’on utilise la seule
observable courbe Force/Déplacement.
Remarque relative au tracé des iso-valeurs : les iso-valeurs de la fonction coût sont tracées
sur la base d'une grille régulière et en utilisant une interpolation linéaire. Les minima locaux
qui apparaissent sur les Figure 4-20 et Figure 4-21 sont le fruit de ce tracé et ne sont pas
représentatifs de multiples extrema.
144
Identification de paramètres matériaux par analyse inverse
fixer ainsi la valeur du paramètre b est une hypothèse forte qui est généralement faite pour
simplifier le problème d'identification. Une telle hypothèse peut cependant fausser les
résultats en manquant un autre jeu de paramètres valable pour une valeur de b différente.
Toute cette problématique vient du fait que l’on essaie d’identifier trois paramètres
d’endommagement sur la base d’une seule observable, globale qui plus est. L’unicité de la
solution n’est alors pas garantie, et il est parfois nécessaire de fixer l’un des paramètres pour
permettre à l’algorithme de minimisation de converger vers un seul minimum. Une autre
solution, plus physique, consisterait à augmenter le nombre d’observables pour enrichir le
problème de minimisation.
3500
3000
2500
2000
F(N)
1500
1000
500 EP
EP + Endommagement
0
0 2 4 6 8 10 12
d(mm)
145
Chapitre 4
Une première identification est menée en ne considérant que deux paramètres à identifier :
ε d et S0. Les autres paramètres du modèle direct sont fixés à leur valeur nominale définie
Tableau 4-2 et rappelée Tableau 4-3. Dans ce même Tableau 4-3, les bornes de recherche des
paramètres à identifier sont définies. L'identification est réalisée grâce à la plateforme MOOPI
en utilisant la méthode de minimisation EGO-p. 10 calculs sont réalisés dans la phase
d'initialisation de la base de données, puis 2 calculs sont réalisés par itération.
L'enrichissement virtuel s'appuie sur la méthode Kriging Believer, décrite au Chapitre 2,
Paragraphe 4.3.
146
Identification de paramètres matériaux par analyse inverse
−1
10
Fc − initialisation
Fc − iteration
Fc − min
−2
Fc
10
−3
10
0 2 4 6 8 10
Iterations
2900
2500
2800
2000
F(N)
F(N)
2700
1500
2600
1000
2500
Observable identifiée Observable identifiée
Observable experimentale Observable experimentale
500 2400
0 2 4 6 8 10 12 4 6 8 10 12
d(mm) d(mm)
(a) (b)
147
Chapitre 4
Surface de Reponse
3
"Vrais" minimum
0.075
Minimum detecté
0.
2.5
07
Iso−fc
5
2 0.0
7 5
5
S0
1.5 0.0 0.05
0.0
1 5
0.
0.05 01
0.5 0.0
0.07 5
5
148
Identification de paramètres matériaux par analyse inverse
« convenables » peuvent parfois être obtenues pour ce type de problèmes d’identification. Cet
extremum supplémentaire, identifié ici, apparaît sur la surface de réponse de la Figure 4-28.
Les paramètres ε d et S0 identifiés (Tableau 4-4) sont bien en accord avec l'analyse par surface
de réponse. Pour une valeur de S0 donnée, la valeur de b a très peu d'influence sur la fonction
coût (Figure 4-16, Figure 4-17), il est donc normal que la valeur de S0 identifiée soit fiable
(erreur inférieure à 1%) malgré l'erreur faite sur b. En revanche les Figure 4-18 et Figure 4-19
montrent une évolution couplée des paramètres ε d et b. Il est donc logique que la valeur
identifiée de ε d soit faussée (erreur de 7.3%), car la valeur de b fixée (b=2.5) ne correspond
pas à celle de la valeur nominale (b=1).
(a) (b)
149
Chapitre 4
Surface de Reponse
1
Minimum detecté
0.9 0.05 Iso−fc
0.05
0.8
0.01
0.7
0.3
0.0 0.05
75
0.2
0.1
0.05 0.1 0.15
Epsd
150
Identification de paramètres matériaux par analyse inverse
et le troisième est fixé à la valeur optimale. Ces surfaces de réponse montrent bien la présence
de minima locaux, dans les plans b - ε d (Figure 4-29a) et ε d -S0 (Figure 4-29c). Ces différents
minima locaux proposent des solutions moins bonnes que la solution nominale, mais
constituent tout de même des solutions acceptables du problème inverse. De plus, le plan de
coupe b - S0 (Figure 4-29b) met à nouveau en évidence la très faible sensibilité de la fonction
coût par rapport au paramètre b autour de la valeur optimale de S0.
Il est important de noter que les surfaces de réponse présentées ici, sont celles construites
durant le processus d'optimisation. Elles ne sont donc pas totalement représentatives des
variations réelles de la fonction, notamment dans les zones où la fonction coût est élevée.
Cette différence de construction, avec les surfaces présentées de la Figure 4-17 à la Figure
4-21, explique les différences observées.
151
Chapitre 4
5
0.0
0.020.03
Minimum detecté Minimum detecté
0.0
0.05
0.04 2
Iso−fc Iso−fc
0.04
0.04
5
0.0
2.5 2.5
0.04
0.03
2 2
0.05
05
0.03
0.
0.05
0.00.504
b
4
0.04
0.0
1.5 4 3 1.5
0.0 0.0
0.0.02
0.
0.0
0
05
3
0.03
1 1 0.
2
1 0.0 05
0.0
0.03
0.02
0.01
0.02
0. 4
03
0.
04
0.5 0.5
0.05 0.1 0.15 0.5 1 1.5 2 2.5 3
Epsd S0
(a) (b)
Surface de Reponse
3
0.
2.5
05
0.04 0.0
2 4
0.03
0.
05
S0
1.5 0.02
0.040.030.02
0.0
0.
03
1 4
0.05
0. 0.0
01
0.04 0.0
0.
0.0 3
02
0.5 5
Minimum detecté
1
Iso−fc
(c)
Figure 4-29 : Surface de réponse - identification virtuelle à trois paramètres - (a) b et ε d - (b) b
et S0 - (c) S0 et εd
152
Identification de paramètres matériaux par analyse inverse
La convergence de la méthode est ralentie (Figure 4-30), par rapport à l'identification des
deux paramètres ε d et S0. Ceci s'explique tout d'abord par le fait que la résolution d'un
problème d'optimisation de dimension plus élevée est plus longue. De plus, l'identification du
paramètre b fait apparaître des minima locaux, ainsi que des zones de sensibilité très faibles,
ce qui a pour effet de ralentir aussi la vitesse de convergence. Après 22 itérations soit 118
calculs de la fonction coût la méthode d'optimisation ne progresse plus.
0
10
Fc − initialisation
Fc − iteration
Fc − min
−1
10
Fc
−2
10
−3
10
0 10 20 30
Iterations
Les tests menés montrent que l'identification des 3 paramètres ( ε d , S0, b) de la loi
d’évolution de l'endommagement est possible malgré la présence de minima locaux et de
zones de sensibilité très faibles. Cependant, tous les tests menés jusqu'à présent faisaient
l'hypothèse que les paramètres de la loi d'écrouissage étaient parfaitement connus. Pour
s'affranchir de cette hypothèse, une identification est menée en considérant les trois
paramètres de la loi d'écrouissage (décrite équation(3.5)) comme des inconnues
supplémentaires.
Afin de limiter la difficulté du problème d'identification, seuls les paramètres ε d et S0 du
modèle d'endommagement sont identifiés. Le paramètre b, qui présente une faible sensibilité,
est fixé à sa valeur nominale. Les plages de variation des différents paramètres sont décrites
dans le Tableau 4-6.
L'identification est réalisée grâce à la plateforme MOOPI en utilisant la méthode de
minimisation EGO-parallèle. 60 calculs sont réalisés dans la phase d'initialisation de la base
de données, puis 8 calculs sont effectués à chaque itération. L'enrichissement virtuel s'appuie
sur la méthode Kriging Believer. L'identification est réalisée sur les 5 paramètres
simultanément.
153
Chapitre 4
Après 240 évaluations du modèle direct, la meilleure solution obtenue est présentée dans
le Tableau 4-6. La fonction coût associée est de 0.93%, La courbe force-déplacement
identifiée est présentée en Figure 4-31.
La solution obtenue montre une erreur inférieure à 5% pour tous les paramètres, excepté
pour la limite élastique. De même que pour le paramètre b, la limite élastique σy semble avoir
peu d’influence sur l’observable. Cela peut provenir également d’une certaine corrélation
entre les paramètres élasto-plastiques, σy et K principalement. Le paramètre σy pouvant être
identifié de manière analytique, il semble préférable de le déterminer au préalable et de le
fixer pour la suite de l'identification.
Remarques :
En plus des performances de l’algorithme de minimisation utilisé, la qualité de
l’identification par analyse inverse dépend de plusieurs facteurs :
o Sensibilité des paramètres sur l’observable : lorsque certains paramètres ont peu
d’influence sur l’observable (comme pour b), leur identification ne peut se faire de
manière précise.
o Corrélation entre paramètres : lorsque deux ou plusieurs paramètres ont des effets
similaires sur l’observable, il est difficile d’identifier un jeu de paramètres unique.
o Nombre de paramètres à identifier : plus le nombre de paramètres à identifier est
important, plus il est difficile de converger de manière précise vers une solution
unique, surtout lorsque l’on n’a qu’une seule observable (globale qui plus est).
Ainsi, l’amélioration du processus d’identification des paramètres nécessite soit de
travailler avec plusieurs observables, soit d’utiliser des observables mieux adaptées
aux paramètres à identifier (peut-être plus locales), pour lesquelles il est plus facile de
découpler les paramètres à identifier.
154
Identification de paramètres matériaux par analyse inverse
Pour tenter de lever cette incertitude, une observable supplémentaire est ajoutée : la courbe
striction/déplacement. L'utilisation de cette observable a été proposée par Mahken [Mah99]
pour identifier les paramètres du modèle d'endommagement de Gurson, dans l'objectif de
lever en partie l'incertitude sur l'unicité de la solution.
On propose ici de mener une analyse similaire à celle proposée pour l'observable force-
déplacement afin de mettre en évidence les apports d'une mesure locale pour l'identification
des paramètres du modèle d'endommagement ductile de Lemaitre.
155
Chapitre 4
Mesure de striction
1.6
1.4
1.2
striction(mm)
0.8
0.6
0.4
0.2
EP
0
0 2 4 6 8 10 12
d(mm)
156
Identification de paramètres matériaux par analyse inverse
(Figure 4-20 pour la fonction coût force-déplacement, Figure 4-39 pour la fonction coût
déplacement/striction).
b
1.4 0.15
1.2 0.06
1.2 0.1
0.04
1
1 0.05 0.02
0.8 0.8
1 2 3 0.4 0.6 0.8 1
S0 S0
0.16
1.6 0.12 5
0.14
0.1
0.12
1.4 4
0.08 0.1
b
0.08
1.2 0.06 3
0.06
0.04
1 2 0.04
0.02
0.02
0.8 1
0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2
S0 Epsd
157
Chapitre 4
S0
1
1.5
0.2
0.8 0.6
1 0.15
0.6
0.4 0.4 0.1
0.5
0.2 0.05
0.2
0 0.1 0.2 0.12 0.14 0.16 0.18 0.2
Epsd Epsd
Une observation plus précise de ces deux zones permet tout de même d'observer une
légère différence sur la forme des surfaces de réponse. Les Figure 4-40 et Figure 4-41
présentent une superposition de la surface de réponse associée à l'observable force-
déplacement et de la surface associée à l'observable déplacement/striction .
La Figure 4-40 montre que l'observable supplémentaire de striction permet de rejeter la
zone S0<0.4 MPa et 0.17< ε d <0.18, zone qui présente une valeur de fonction coût force-
déplacement très faible. La fonction coût déplacement/striction n'apporte pas d'information
supplémentaire dans le reste du domaine de recherche, sur ce plan de coupe.
1.2
Déplacement/Force
Déplacement/Striction
1
Minimum global
0.8
S0
0.4
0.2
0.12 0.14 0.16 0.18 0.2
Epsd
158
Identification de paramètres matériaux par analyse inverse
La Figure 4-41 détaille la forme des deux surfaces de réponse dans la zone du second
minimum détecté dans le plan b- ε d (le minimum local non désiré). Les iso-valeurs de la
fonction coût striction-déplacement (Figure 4-41, trait plein) présentent une zone minimum
décalée par rapport aux iso-valeurs de la fonction coût force-déplacement (Figure 4-41, trait
discontinu). La combinaison de ces deux fonctions coût permet, d'après l'analyse menée ici,
d'atténuer l'effet des minima locaux qui apparaissent ici non coïncidents.
3
Déplacement/Force
Déplacement/Striction
2.8
2.6
b
2.4
2.2
2
0.1 0.12 0.14 0.16
Epsd
Il faut cependant bien noter que l'analyse de sensibilité menée ne permet pas de conclure
que la mesure de la striction permet de lever le problème de l'unicité de la solution. Nous
retrouvons toujours une faible sensibilité par rapport au paramètre b.
L'identification est réalisée dans les mêmes conditions que celles décrites Tableau 4-5. Les
résultats obtenus sont présentés dans le Tableau 4-7. La valeur de la fonction coût minimum
est de 0.22%, avec respectivement fc force = 0.17% et fcstriction = 0.43% .
159
Chapitre 4
Cependant, l'ajout de l'observable de striction permet tout de même de limiter les minima
locaux. Cet effet est notamment visible sur la Figure 4-43.b, où le minimum local présent en
b=2.5 sur la surface de réponse associée à la fonction coût force-déplacement (Figure 4-29 b),
disparaît.
Cet ajout de l'observable de striction permet également une convergence plus rapide de la
méthode de minimisation. Une solution de qualité équivalente en terme de valeur de fonction
coût sur l'observable en force est trouvée en 10 itérations de la méthode, contre 22 pour
l'identification réalisée exclusivement avec l'observable en force (Figure 4-30).
0
10
Fc − initialisation
Fc − iteration
Fc − min
−1
10
Fc
−2
10
−3
10
0 2 4 6 8 10
Iterations
Figure 4-42 : Convergence - identification virtuelle à trois paramètres - fonctions coût force-
déplacement et striction-déplacement
160
Identification de paramètres matériaux par analyse inverse
0 .0 6
0.0 Minimum detecté
0.090.06
3
0.09
Iso−fc
0.06
2.5 2.5
2
0.06
09
0.1
0.12
0.
2 2
0.03
0.06
0.06
0.03
0.03
0.09 .06
b
0.03
0.1
0
0.09
1.5 1.5 2
0.06
6 0.00
0.0 6
0.12
1 1
0.12
0.03
0.
0.06
0.06
09
0.03
0.03
Minimum detecté
0.09
0.1
0.0
Iso−fc
0.0
0.
2
9
6
03
0.5 0.5
0.05 0.1 0.15 0.5 1 1.5 2 2.5
Epsd S0
(a) (b)
Surface de Reponse
3
0.09 2
0.1
2.5
0.12
2 0.09
0.06 0.06 0.
0.03 09
S0
1.5
0.0
0.06
0.03
6
0.06
1 0.09
0.0
0.03
3
0.1 0.09
0.15 2
0.5 0.
0.1 12detecté
Minimum 0.06
Iso−fc 5
(c)
Figure 4-43 : Surfaces de réponse - identification virtuelle à trois paramètres - fonctions coût
force-déplacement et striction-déplacement - (a) b et ε d - (b) b et S0 - (c) ε d et S0
161
Chapitre 4
converge vers la valeur nominale de 0.16. En revanche pour le paramètre b il n'y a pas de
convergence vers une solution unique. La Figure 4-44c fait apparaître 3 zones de convergence
différentes. L'une de ces zones correspond bien à la valeur nominale du paramètre b=1, mais
d'autres valeurs de b semblent donner des solutions tout aussi satisfaisantes (b=0.5). Ces
graphiques mettent là encore en évidence la difficulté d'identifier le paramètre b du modèle
d'endommagement de Lemaitre.
3.5 0.2
S0 EpsD
0.18
3
0.16
2.5 0.14 Convergence
2 0.12
0.1
1.5
0.08
Convergence
1 0.06
0.04
0.5
0.02
0 Fc 0 Fc
0 0.02 0.04 0.06 0.08 0 0.02 0.04 0.06 0.08
(a) (b)
3.5
b
3
Minimum Local
2.5
Minimum Local
0.5
0 Fc
0 0.02 0.04 0.06 0.08
(c)
162
Identification de paramètres matériaux par analyse inverse
problèmes liés à la sensibilité de certains paramètres restent présents. On réalise tout de même
une identification virtuelle des paramètres d'endommagement et de plasticité. Cette
identification est réalisée dans les mêmes conditions que celles présentées au paragraphe
3.3.5. Ces conditions, ainsi que les résultats sont présentés Tableau 4-8.
Apres 130 évaluations du modèle direct, la fonction coût obtenue est de 1.07%. Les deux
observables sont présentées Figure 4-45et Figure 4-46.
Les résultats obtenus peuvent être considérés comme bons d'après la valeur de la fonction
coût et l'allure des observables identifiées. Les paramètres identifiés sont cependant identifiés
avec une erreur relativement importante, en particulier sur le paramètre σy. L'ajout de
l'observable de striction dans le calcul de la fonction coût ne permet donc pas une
identification plus précise, mais accélère tout de même énormément la vitesse de convergence
de la méthode de minimisation. En effet la même identification sans prendre en compte
l'observable de striction propose une solution convenable après 240 évaluations du modèle
direct contre 130 pour l'identification tenant compte de la striction.
Cette accélération peut nous laisser penser que l'ajout de l'observable de striction permet
de mieux conditionner le problème de minimisation et de limiter le nombre des minima
163
Chapitre 4
164
Identification de paramètres matériaux par analyse inverse
F
Obsexp
Obsnum
dmax d
(a) phase 1
F
Obsexp
Obsnum
(b) phase 2
165
Chapitre 4
166
Identification de paramètres matériaux par analyse inverse
Nombre de
Plage de Paramètres Fonction
Phase Paramètres calculs du modèle
recherche identifiés coût
direct
σ y (MPa) [10 400] 93 MPa 139 calculs du
1 K (MPa) [200 1000] 441 MPa 0.77% modèle direct /
27 itérations
n [0.1 0.6] 0.44
b [0.5 3] 1.00 104 calculs du
2 S0 [0.1 3] 0.77 MPa 0.77% modèle direct /
εd [0.04 0.20] 0.097 18 itérations
Tableau 4-9 : Bornes de recherche et résultat de l'identification des paramètres de l'aluminium
[AlMg3]
167
Chapitre 4
168
Identification de paramètres matériaux par analyse inverse
De nouveaux modèles prenant en compte ces deux paramètres [Bai08] ainsi que de
nouvelles procédures d’identification ont été proposés. Les procédures d’identification
montrent qu’il est important d’identifier les paramètres du modèle d’endommagement dans
des conditions (en terme de trajet de chargement) les plus proches possible du procédé de
mise en forme étudié.
Ce n’est pas le cas ici entre notre essai de traction et la sollicitation subie par la matière
dans la zone endommagée lors du clinchage. La Figure 4-51 illustre le trajet de chargement
que subit la matière pendant le procédé de clinchage et celui subi par la matière pendant
l'essai de traction. Ce graphique, tracé dans l’espace Triaxialité des contraintes Tx / Angle de
Lode θ , montre clairement que l'identification des paramètres matériau en traction n'est pas
réalisée dans les mêmes conditions que le procédé de clinchage.
Cette problématique est traitée actuellement au CEMEF dans le cadre de la thèse de J.M.
Gachet [Gac10].
Tx
Figure 4-51 : Représentation de l'état des contraintes dans le plan de la triaxialité des
contraintes et de l’angle de Lode [Bai08] - trajet de chargement du procédé de clinchage et de
l'essai de traction
169
Chapitre 4
4 Conclusion
Dans ce chapitre, les principes généraux de l'identification de paramètres matériau par
analyse inverse ont été présentés dans un premier temps. Nous avons montré comment
l'analyse inverse s'inscrit dans le cadre, plus général, de la problématique d'optimisation. Cette
mise en relation, analyse inverse / optimisation, permet de s'appuyer sur la plateforme
MOOPI pour réaliser l'identification.
Un premier pas a été réalisé pour tenter de lever ces incertitudes en investiguant, par une
approche numérique, l'intérêt de l'ajout d'une observable supplémentaire : une mesure de
striction dans la section utile de l'éprouvette. Les résultats obtenus sont en demi-teinte. En
effet, l'ajout de l'observable de striction ne procure pas vraiment une meilleure précision de la
solution identifiée. Les problèmes de faible sensibilité des observables vis-à-vis de certains
paramètres subsistent. Il est à noter également que cette mesure de striction supplémentaire
permet d’enrichir plutôt l'identification de la partie plasticité que de la partie
endommagement, ce qui explique en partie la faible amélioration de l’identification du
paramètre d’endommagement b en rajoutant cette observable locale. Cependant, l'ajout de
cette observable permet de limiter les minima locaux dans le domaine de recherche et donc
d'accélérer de manière significative la vitesse de convergence de la méthode de minimisation.
L'ajout d'une information locale présente donc un intérêt.
170
Identification de paramètres matériaux par analyse inverse
sûr plus sensible aux paramètres d’endommagement que le début de cette courbe.
L’adaptation du calcul de la fonction coût à la rupture en fin d’essai est aussi un
paramètre important à prendre en compte.
• Représentativité des lois de comportement et du modèle d’endommagement :
quelle que soit la méthode de minimisation choisie, ou le nombre d’observables
utilisé, la procédure d’identification ne peut pas fournir une représentation correcte
de l’expérience si le modèle de comportement, ou le modèle d’endommagement
n’est pas adapté au matériau étudié.
• Corrélation entre certains paramètres des lois : comme nous l’avons vu, certains
paramètres du modèle ont des effets corrélés sur l’observable. Il est alors difficile
de les identifier de manière unique. Il faut dans ce cas les identifier de manière
séparée, à partir d’observables plus sensibles à l’un ou à l’autre de ces paramètres.
• Problème d'unicité des paramètres d'endommagement : des études récentes ont
montré la sensibilité de l'endommagement ductile vis-à-vis de la triaxialité des
contraintes Tx et de l'angle de Lode θ [Bai08]. Une autre manière de converger
vers une solution unique est de baser l'identification sur plusieurs essais
mécaniques faisant varier Tx et θ . Cette approche n'a pas été traitée ici mais est
abordée dans le cadre du travail de thèse de J.M. Gachet [Gac10].
171
Identification : apport des mesures de champs
Chapitre 5
Identification : apport des mesures de champs
173
Chapitre 5
174
Identification : apport des mesures de champs
Problème
physique
Modèle EF
Essais
(
FEF P, Uˆ )
Fonction coût
( )
Mesures de champs
R f ( P ) = FEF P, Uˆ − Fˆ Accord
F̂ et Û
Solution du problème
d'identification
Figure 5-1 : Schéma de résolution de la méthode par recalage éléments finis REF-F
L’autre méthode, REF-U, est une formulation en déplacement. Cette méthode alternative,
basée sur la concordance des champs de déplacement, consiste à minimiser le système
suivant :
(
P = arg min U ( P ) − Uˆ ) (5.2)
Où U(P) est le champ obtenu par le calcul éléments finis.
La minimisation se fait alors sous contrainte d'égalité entre la force globale mesurée et
celle obtenue par le calcul éléments finis.
Il existe ensuite plusieurs approches généralisées mêlant les méthodes REF-F et REF-U.
Giton et al. [Git06] et Pagnacco et al. [Pag07] proposent des formulations généralisées, par
exemple la formulation multi-objectifs suivante :
(
P = arg min R ( P ) = U ( P ) − Uˆ , R ( P ) = F (U , P ) − Fˆ
U F int ) (5.3)
En conclusion, l’identification par recalage éléments finis permet, d’une part de traiter des
géométries complexes, et d’autre part l’utilisation d’un champ de déplacements mesurés (il
n’est pas nécessaire de connaître le champ de déformations). Néanmoins, la mise en œuvre
d’une telle méthode nécessite de choisir une formulation du problème inverse ainsi qu’une
modélisation pour les conditions aux limites.
Concernant les méthodes présentées ci-dessus, la formulation en effort présente l’avantage
de ne pas être sensible aux mouvements de corps rigides, mais elle nécessite une parfaite
connaissance du champ de déplacement et elle s’avère assez sensible aux bruits de mesure. La
formulation en déplacement est quant à elle moins sensible aux bruits de mesure. Elle permet
d’utiliser des mesures incomplètes en déplacement, mais en contre partie elle s’avère être
sensible aux mouvements de corps rigides et elle est toujours résolue de manière itérative.
175
Chapitre 5
La formulation généralisée est plus souple. Les mesures cinématiques ne sont utilisées
qu’aux bords. Cette approche est donc particulièrement intéressante en 3D avec des mesures
de surface.
176
Identification : apport des mesures de champs
Nous ne développerons pas davantage cette méthode, qui est plus une méthode servant à
détecter la présence de singularités qu'une méthode permettant d'identifier des paramètres
matériaux.
Cette expression est vraie pour tout champ virtuel cinématiquement admissible. De plus, il
existe une équation pour chaque champ virtuel. Pour un problème donné, les équations
disponibles sont les suivantes :
• Equation d’équilibre :
Forme forte :
σ ij , j + fi = 0 + conditions aux limites (5.6)
Ou l’expression écrite en formulation faible :
−Cijkl ∫ ε kl ε ij* dV + ∫ Ti ui*dS + ∫ fi ui*dV = 0 (5.7)
V S V
σ ij = Cijkl ε kl (5.8)
• Lien déplacements-déformations
1
ε ij =
2
( ui, j + u j ,i ) (5.9)
A partir des équations précédentes on peut écrire un système de P équations à P
inconnues :
−Cijkl ε kl ε ij(1)*dV + Ti ui(1)*dS = 0
∫
V
∫S
−Cijkl ∫ ε kl ε ij dV + ∫ Ti ui dS = 0
(2)* (2)*
V S (5.10)
...
−Cijkl ∫ ε kl ε ij( p )*dV + ∫ Ti ui( p )*dS = 0
V S
177
Chapitre 5
178
Identification : apport des mesures de champs
Le champ de déplacement ainsi obtenu est disponible sur une face, avec une résolution
temporelle dépendante de la fréquence d'échantillonnage de la caméra, et une résolution
spatiale fixée par la taille des fenêtres de corrélation et leur taux de recouvrement. La
précision de la mesure est liée notamment à la qualité du mouchetis réalisé sur la face de
l'éprouvette. Le mouchetis de peinture doit être le plus fin possible et proposer le plus large
spectre de nuance de gris possible afin d'obtenir une corrélation de qualité.
Essai mécanique
Corrélation d'image
Champs de déplacement
Mouchetis
Les champs de déplacement ainsi obtenus sont convertis dans un format lisible par la
librairie CimLib de calcul éléments finis. Cette conversion est réalisée via une série de scripts
Matlab®. Trois grandeurs sont alors disponibles sur une grille régulière (correspondant aux
fenêtres de corrélation utilisées) : le déplacement selon X, selon Y et un champ binaire sur la
disponibilité de la mesure (Figure 5-3). Ces données sont disponibles pour les différents
instants fixés par la fréquence d'acquisition des images.
179
Chapitre 5
La grille de support de ces mesures est maillée par des éléments triangles pour pouvoir
utiliser les opérateurs de transport disponibles dans la librairie CimLib [Dig07] qui travaille
exclusivement avec des simplexes (triangles, tétraèdres).
La méthode mise en place ici impose le pas de temps du calcul par éléments finis. En
effet, celui-ci est imposé par la fréquence d'acquisition de la prise de vue. Cette procédure est
donc itérative, elle est schématisée Figure 5-5.
180
Identification : apport des mesures de champs
Figure 5-4 : Principe du pilotage d’une simulation éléments finis avec des mesures de champs
Calcul mécanique EF
ti+1 = ti+dt
Figure 5-5 : Processus de pilotage d'un calcul éléments finis par les mesures de champs
181
Chapitre 5
2.2.2 Validations
La superposition des champs issus de mesures de champs (Figure 5-6, zone large) et ceux
issus du calcul éléments finis (Figure 5-6, zone centrale) permettent de valider la
méthodologie proposée. Les iso-valeurs de déplacement se superposent bien sur le bord où les
conditions limites sont imposées (Figure 5-6).
(a) iso déplacement suivant x (t=23s) (b) iso déplacement suivant y (t=23s)
Calcul EF Calcul EF
(c) iso déplacement suivant x (t=114s) (d) iso déplacement suivant y (t=114s)
Figure 5-6 : Validation conditions limites issues des mesures de champs - superposition des
iso-valeurs de déplacement issues des mesures de champs et obtenues par calcul éléments finis
182
Identification : apport des mesures de champs
2.3.1 Formulation
L'étape suivante dans l'exploitation des mesures de champs est de mettre en place une
fonction coût adaptée aux mesures de champs, c'est-à-dire qui est capable d'évaluer l'écart
entre une carte de déplacement numérique et expérimentale.
Pour réaliser cette évaluation on dispose, à chaque incrément, d'une carte de déplacement
expérimentale rattachée au maillage régulier issu de la prise de vue, et d'une carte de
déplacement numérique rattachée au maillage courant du calcul éléments finis.
Il est primordial de comparer les déplacements numériques et expérimentaux du même
point matériel. La première solution envisagée est de travailler avec un maillage coïncidant
entre les points de mesure expérimentaux et les nœuds du maillage éléments finis. Cette
approche est proposée par Springmann et al. [Spr06]. L'écart entre les déplacements
numériques et expérimentaux est alors évalué directement nœud à nœud. Un exemple de
maillage utilisé par Springmann est présenté Figure 5-7.
Figure 5-7 : Maillage coïncident entre les points de mesure du déplacement et les nœuds du
maillage éléments finis [Spr06]
Cette approche est peu flexible. En effet, le maillage est directement fixé par les données
expérimentales. Pour s'affranchir de cette contrainte et afin de travailler avec un maillage
quelconque, il est nécessaire de comparer les champs de déplacement rapportés à la
configuration initiale. La démarche mise en place est décrite Figure 5-8.
183
Chapitre 5
Uimp
Uimp
Calcul EF
Unum
U −1
Unum
Transport
Uimp
Uimp
Uimp Uimp
U −1
Essai - MDC Uexp Uexp Uexp, Unum
184
Identification : apport des mesures de champs
fcMDC ( P ) = ω x fcx ( P ) + ω y fc y ( P )
avec
∑ δ . (Ux ( P ) − Ux )
n 2
T
num mdc
i i i
fcx ( P ) = ∑ i =1
n
t =1
∑δ i
i =1
n 2
T ∑ (
δ i . Uyinum ( P ) − Uyimdc
)
fc y ( P ) = ∑ i =1
n
(5.11)
t =1
∑ δi
i =1
où T est le nombre de clichés, n est le nombre de points de mesure, δ i la valeur binaire sur
la disponibilité de la mesure valant 0 (mesures non disponibles) ou 1(mesures disponibles),
ω x et ω y la pondération donnée aux composantes du déplacement, Uxinum et Uyinum les
déplacements suivant les deux axes issus du calcul éléments finis, et Uximdc et Uyimdc ceux issus
des mesures de champs. Afin de ne pas surcharger les notations, l'indice d'incrémentation
temporel t n'a pas été écrit : les grandeurs Uxinum , Uyinum , Uximdc , Uyimdc et δ i sont dépendantes
de t.
La formulation proposée ici permet donc de comparer et d'exploiter des mesures de
champs dans le cadre de grandes déformations et en utilisant toutes les méthodes numériques
spécifiques à ce cas, notamment les outils de remaillage adaptatif. Cette formulation ne tient
pas compte des incertitudes de mesure. Les travaux de Fazini [Faz09] proposent une
formulation proche de celle décrite équation (5.11), qui permet d'intégrer une incertitude
constante ou proportionnelle dans l'évaluation de la fonction coût.
La formulation (5.11) peut être rattachée aux méthodes de recalage éléments finis en
déplacement (REF-U) présentées précédemment. Cependant il semble préférable de conserver
l'information de la force de traction, la fonction coût utilisée sera donc une somme pondérée
des fonctions coût liées à l'observable globale de force et à l'observable locale de mesure de
champs de déplacement.
185
Chapitre 5
Uxinum = Uxinum
,t −1
Si δ inum =
,t
,t 0, alors num
Uyi ,t = Uyi ,t −1
num
(5.12)
,t = Uxi ,t −1
Uxiexp exp
Si δ = 0, alors exp
exp
i ,t
Uyi ,t = Uyi ,t −1
exp
35
Une simulation d'un essai de traction sur cette éprouvette est réalisée, les informations
(force résultante et champs de déplacement) obtenues sont conservées. Elles serviront de
« données expérimentales virtuelles » pour l'analyse de sensibilité. Les paramètres matériaux
utilisés pour la génération de l'essai virtuel sont rappelés dans le Tableau 4-2.
186
Identification : apport des mesures de champs
Paramètres Jeu 1
E 69000 MPa
ν 0.3
σy 46 MPa
K 430MPa
n 0.34
b 1
S0 0.7 MPa
wc 0.8
εd 0.16
lc 0.06 mm
h 0.2
Tableau 5-1 : Jeu de paramètres de référence pour la génération de l'essai virtuel
Le modèle direct mis en place ne modélise que la partie centrale de l'éprouvette (une zone
de 10 mm de part et d'autre du plan moyen de l'éprouvette), comme illustré Figure 5-4. On
utilise ici les mesures de champs à double titre : pilotage en déplacement sur les bords de la
zone d’étude et évaluation de la fonction coût (équation(5.11) ).
187
Chapitre 5
Figure 5-10: Surfaces de réponse - analyse de sensibilité des contributions de la fonction coût par rapport aux paramètres de la loi d'écrouissage
188
Identification : apport des mesures de champs
Tout d'abord, nous pouvons noter que la forme des surfaces de réponse associées à la
contribution en force (Figure 5-10, colonne 1) ont une allure classique : le problème est mal
conditionné, les surfaces de réponse présentent une vallée importante. Ces vallées sont
significatives d'un problème de corrélation entre les paramètres du modèle, notamment
lorsque l’on n’utilise qu’une observable globale.
Plusieurs remarques peuvent être formulées par rapport à l'apport des mesures de champs :
• Analyse de sensibilité par rapport aux couples de paramètres (K, n) et (K, Sig0)
(Figure 5-10, lignes 1 et 2) :
o Observable force (Figure 5-10, colonne 1, lignes 1 et 2) : les deux surfaces
de réponse présentent une vallée relativement étroite (pointillés rouges), un
ensemble de jeux de paramètres (K, n) et (K, Sig0) proposent des fonctions
coût de valeur faible et sont donc des solutions potentielles du problème.
L'identification menée sur cette seule observable présente donc une
incertitude importante, due à la corrélation entre les paramètres du modèle.
o Observable Ux (Figure 5-10, colonne 2, lignes 1 et 2) : les deux surfaces de
réponse présentent une zone de minimum relativement large. La sensibilité
de la fonction coût aux paramètres est donc relativement faible. On
remarque cependant que les vallées observées précédemment ont disparues.
Ceci dénote d’une meilleure décorrélation des paramètres.
o Observable Uy (Figure 5-10, colonne 3, lignes 1 et 2) : les deux surfaces de
réponse présentent des vallées relativement étroites (pointillés rouges). Ces
vallées dénotent encore une fois d’une certaine corrélation entre les
paramètres vis-à-vis de l’observable. Cependant, ces vallées ont comme
principale caractéristique d'être orientées différemment de celles observées
sur la surface de réponse relative à la force (Figure 5-10, colonne 1, lignes
1 et 2) et sont donc, à ce titre, très riches d'informations.
o Combinaison des 3 contributions (Figure 5-10, colonne 4, lignes 1 et 2) :
une combinaison de ces surfaces de réponse permet de lever les problèmes
de sensibilité et de corrélation des paramètres liés aux couples de
paramètres (K, n) et (K, Sig0). En effet, ces surfaces de réponse présentent
un minimum bien défini et ne font plus apparaitre de vallées,
caractéristiques d'une corrélation entre les paramètres du modèle.
• Analyse de sensibilité par rapport au couple de paramètres (n, Sig0) (Figure 5-10,
ligne 3) : les trois surfaces de réponse (Figure 5-10, ligne 3, colonnes 1, 2 et 3)
présentent des vallées similaires (pointillés rouges). Pour ce couple de paramètres
les mesures de champs ne permettent pas de lever les problèmes de corrélation
entre les paramètres n et Sig0. Ce problème de corrélation est récurrent pour les
différentes observables, l'une des solutions couramment employée est de fixer le
seuil Sig0 en analysant la courbe force déplacement.
Cette analyse de sensibilité par rapport aux paramètres de la loi d'écrouissage, réalisée sur
un essai virtuel, permet de mettre en évidence l'intérêt d'exploiter des observables issues de
mesures de champs : cet ajout permet d'obtenir un problème de minimisation mieux
conditionné, limitant les corrélations entre paramètres. Le problème d'inversion est donc
mieux posé, il sera donc plus simple à résoudre et les paramètres identifiés le seront avec plus
de fiabilité.
189
Chapitre 5
190
Chapitre 5
Figure 5-11: Surfaces de réponse - analyse de sensibilité des contributions de la fonction coût par rapport aux paramètres de la loi d'endommagement
192
Identification : apport des mesures de champs
• Analyse de sensibilité par rapport au couple de paramètres (b, EpsD) (Figure 5-11,
ligne 3) :
o Observable force (Figure 5-11, ligne 3, colonne 1) : la surface de réponse
présente une vallée (pointillés rouges), tout de même moins prononcée que
sur les surfaces préalablement étudiées.
o Observables Ux et Uy (Figure 5-11, ligne 3, colonnes 2 et 3) : les surfaces
de réponse présentent également des vallées, qui sont cependant moins
prononcées que pour l’observable Force.
o La somme des trois surfaces de réponse (Figure 5-11, ligne 3, colonne 4)
fait toujours apparaitre cette vallée. Elle est cependant un peu moins
prononcée que celle relative à l'observable de force. Le paramètre EpsD
semble clairement identifié, mais c’est le paramètre b qui présente une zone
(0.5<b<1.1) de faible sensibilité de l’observable par rapport à ce paramètre.
Cette analyse de sensibilité permet de montrer l'apport des mesures de champs pour
l'identification des paramètres de la loi d'évolution de l'endommagement.
La faible sensibilité par rapport au paramètre b (à S0 fixé) reste présente mais moins
marquée.
193
Chapitre 5
L’amélioration obtenue ici est en accord avec l’amélioration qui avait été obtenue avec la
mesure locale de la striction. Sur le plus long terme, la mesure des champs de déplacement
peut s’avérer plus intéressante, notamment pour l’identification de propriétés anisotropes.
Cette analyse de sensibilité met aussi en évidence une difficulté liée aux mesures de
champs et plus particulièrement à l'aspect multi-objectifs de l'identification. La fonction
minimisée est une somme pondérée de 3 contributions. Or pour le cas traité ici les 3
contributions ont un ordre de grandeur différent (une décade entre la contribution de la force
et les contributions des mesures de champs). La pondération choisie aura donc un impact fort
sur la solution du problème.
194
Identification : apport des mesures de champs
Pour compléter ces mesures, on propose de réaliser une extrapolation des mesures suivant
des lignes perpendiculaires à l'axe de l'éprouvette (Figure 5-12 a et b, lignes blanches). Les
mesures disponibles sur ces lignes (Figure 5-13, cercles bleus) sont ensuite extrapolées par un
polynôme du second ordre, polynôme qui sert ensuite à compléter les mesures (Figure 5-13,
points rouges). Cette méthode de complétion des mesures est appliquée dans les zones
supérieure et inférieure de l'éprouvette, zones qui serviront ensuite pour le pilotage du modèle
éléments finis. Les cartes de déplacement ainsi obtenues sont présentées Figure 5-14. Il est à
noter que ces mesures complétées ne servent que pour le pilotage de la simulation ; seules les
données brutes sont utilisées pour le calcul de la fonction coût.
195
Chapitre 5
Figure 5-13 : Extrapolation des mesures sur le bord de l'éprouvette pour l’application des
conditions aux limites
196
Identification : apport des mesures de champs
On dispose, pour cet essai, de 140 clichés allant du début de l'essai à la rupture complète
de l'éprouvette, ainsi que de l'enregistrement de l'effort de traction (Figure 5-15).
Figure 5-15 : Effort de traction sur éprouvette de traction DC04 - Identification de la loi
d'écrouissage
197
Chapitre 5
• Pour ce test de faisabilité un seul essai a été réalisé. La fiabilité des mesures
expérimentales n'est donc pas assurée. Des essais de répétabilité sont donc
indispensables.
(a) Ecart sur le champ de déplacement selon x (b) superposition des champs de
- t=84s déplacement suivant x expérimentaux (trait)
et numériques (couleur) - t=84s
(c) Ecart sur le champ de déplacement selon y (d) superposition des champs de
- t=84s déplacement suivant y expérimentaux (trait)
et numériques (couleur) - t=84s
198
Identification : apport des mesures de champs
199
Chapitre 5
Figure 5-17 : Identification de la loi de plasticité de l'acier DC04 - surfaces de réponse - contributions des fonctions coût
200
Conclusion générale et perspectives
Le tracé des surfaces de réponse présentées sur la Figure 5-17 permet de faire plusieurs
remarques :
• Surfaces de réponse relatives à la force (Figure 5-17, colonne 1) : la forme de ces
trois surfaces est comparable à celle obtenue lors de l'étude de sensibilité. Des
vallées, caractéristiques de la corrélation des paramètres vis-à-vis de l’observable
force, sont visibles. Sur chacune de ces surfaces, le minimum détecté par rapport à
la fonction coût globale (donc le résultat final de toute la procédure
d’identification) est représenté par une croix bleue.
Les résultats obtenus sont présentés sur les Figure 5-18 et Figure 5-19, et les paramètres
identifiés sont donnés dans le Tableau 5-3.
201
Conclusion générale et perspectives
Figure 5-18 : Effort de traction sur éprouvette de traction DC04 - Identification de la loi
d'endommagement
202
Conclusion générale et perspectives
(a) Ecarts sur le champ de déplacement (b) Ecarts sur le champ de déplacement
selon y - t=153s selon x- t=153s
Les résultats proposés ici sur l'identification des paramètres de la loi d'endommagement à
partir des mesures de champs ne sont pas vraiment satisfaisants. Le fort écart déjà introduit
en termes de déplacements lors de l'identification des paramètres du modèle d'écrouissage
peut expliquer l'inefficacité de l'identification des paramètres de la loi d'endommagement sur
cette base.
5 Conclusion
Dans ce chapitre, une approche permettant d’exploiter les mesures de champs pour
l'identification des paramètres matériaux a été mise en place dans le cadre de grandes
203
Conclusion générale et perspectives
déformations. L'approche retenue pour notre étude est une approche par recalage éléments
finis. La méthodologie appliquée est basée sur deux étapes : le pilotage du modèle direct par
des mesures issues des mesures de champs, et la construction d'une fonction coût comparant
les champs de déplacement numériques et expérimentaux.
L'analyse de sensibilité par surface de réponse, menée sur un essai virtuel, montre l'apport
très intéressant des mesures de champs pour identifier les paramètres de loi d'écrouissage : la
plupart des problèmes de corrélation entre les paramètres sont levés. Les gains sont moins
importants sur l'identification des paramètres de la loi d'évolution de l'endommagement, mais
on note tout de même que le couple de paramètres (S0, EpsD) est plus facilement identifiable.
La faible sensibilité par rapport au paramètre b reste présente.
204
Conclusion générale et perspectives
Depuis quelques années, il a été démontré l’intérêt d’intégrer la phase de mise en forme
dans le dimensionnement des structures mécaniques en service. Ainsi, l’optimisation de la
tenue mécanique des structures nécessite la construction d’une chaîne de simulation virtuelle
intégrant mise en forme et propriétés en service. C’est dans le cadre de l’assemblage par
déformation plastique que cette approche a été illustrée ici.
Ce travail de thèse a permis dans une première partie de mettre en place la plateforme
MOOPI (MOdular software dedicated to Optimization and Parameters Identification). Cet
outil, dédié à l'étude de sensibilité, l'analyse inverse et l'optimisation, est capable de
communiquer avec un large panel de logiciels de calcul scientifique. La solution proposée
atteint les objectifs initialement fixés : la plateforme est capable de résoudre différents types
de problèmes d'optimisation et sa structure orientée objet permet l'intégration de nouvelles
méthodes d'analyse et d'optimisation, ce qui assure son évolutivité.
Un algorithme d'optimisation adapté à la problématique des assemblages par déformation
plastique a été implémenté. Une version parallèle de l'algorithme est également proposée. Les
principales caractéristiques et points forts de l'algorithme EGO-p implémentés dans MOOPI
sont les suivants :
• Une méthode d’optimisation globale, évitant de rester piégée dans des
minima locaux.
• Un algorithme adapté aux problèmes lourds en terme de temps de calcul
(optimisation assistée par métamodèle par krigeage et parallélisation de la
méthode).
• Un faible nombre de paramètres de réglage de l'algorithme, ce qui facilite son
utilisation dans un cadre industriel.
• Un algorithme non intrusif par rapport au logiciel de calcul, permettant ainsi de
coupler MOOPI à différents logiciels de calcul.
La seconde partie du travail est axée sur l'exploitation de la plateforme MOOPI. Deux
applications principales ont été développées :
• La mise en place d'une procédure d'optimisation des procédés d'assemblage par
déformation plastique.
• L'identification par analyse inverse des paramètres d'une loi de comportement
élasto-plastique couplée à un modèle d’endommagement ductile de Lemaitre.
205
Conclusion générale et perspectives
La seconde application fait l'objet des chapitres 4 et 5. Dans ces chapitres ont été
investiguées les possibilités d'identification des paramètres de la loi d'endommagement
ductile de Lemaitre sur la base d'essais de traction. La définition d'une fonction coût adaptée
au problème d'adoucissement/rupture a permis tout d'abord de mieux définir le problème
inverse. Ensuite, une analyse de sensibilité par étapes successives a été réalisée : à chaque
étape la base d'observable a été enrichie. Tout d'abord l'analyse de sensibilité menée par
rapport à l'observable force/déplacement a permis de mettre en évidence la difficulté liée à
l'identification du paramètre b du modèle d'endommagement de Lemaitre. En effet les
surfaces de réponse tracées font apparaitre une très faible sensibilité de la fonction coût par
rapport à ce paramètre ainsi que des minima multiples. Ensuite, la base d'observables a été
enrichie par une mesure locale de striction de l'éprouvette. Contrairement à ce qui était
attendu, cet ajout ne modifie pas sensiblement la nature du problème inverse : les faibles
sensibilités et les minima locaux restent présents. Cependant les sensibilités sont améliorées,
ce qui explique l'accélération de la vitesse de convergence de la minimisation. Cette
accélération met en évidence l'intérêt d'enrichir la base d'observables.
Finalement, au chapitre 5, des mesures de champs sont ajoutées à la base d'observables.
Les tests numériques menés montrent que ce type d'observables est très riche et améliore la
forme du problème inverse. Cette amélioration est très nette pour l'identification des
paramètres de la loi d'écrouissage : les problèmes de corrélation entre paramètres sont levés
grâce à l'utilisation des mesures de champs. Concernant les paramètres d’endommagement, la
corrélation entre les paramètres S0 et ε d est aussi diminuée. En revanche, même si l’on peut
noter une légère amélioration, l'identification du paramètre b reste délicate. On remarque
notamment l’importance du choix des poids associés aux différentes parties de la fonction
coût. Le premier test d'identification à partir de mesures de champs expérimentales n'est pas
tout à fait concluant, et plusieurs raisons ont été évoquées : précision des mesures de champs
expérimentales, représentativité de la loi de comportement élasto-plastique utilisée,
anisotropie… Ce test a cependant permis de tester la totalité de la méthodologie qui a été mise
en place, notamment le traitement du bruit de mesure et le lien entre le logiciel de corrélation
d'image (Aramis®) et la librairie CimLib®.
Au travers de ces deux derniers chapitres, il est important de noter la grande difficulté de
dissocier la plasticité de l’endommagement ductile à partir d’observables globales (courbes
force-déplacement) et locales (mesures de striction). L’observation des mécanismes
d’endommagement à l’échelle de la microstructure peut permettre, à terme, de mieux
206
Conclusion générale et perspectives
Globalement, ces travaux ont permis de mettre en place un outil d'optimisation efficace et
fiable, qui a pour vocation d'être utilisé au sein du laboratoire pour d’autres applications.
Plusieurs collaborations ont d’ailleurs été mises en place au cours de la thèse pour utiliser la
plateforme MOOPI avec d’autres logiciels ou pour d’autres problèmes physiques. Quelques
exemples sont présentés en Annexe pour illustrer cela. Dans le cadre du projet MONA LISA,
cet outil a permis d'apporter des solutions en termes d'identification de paramètres matériaux
et d'optimisation des assemblages par déformation pastique. Cependant de nombreuses
perspectives et améliorations sont envisageables.
207
Conclusion générale et perspectives
mesurer des déplacements dans des conditions extrêmes, aussi bien en grande déformation
qu’à haute température. La Figure i, montre l’évolution du champ de déplacement, et du
champ de déformation sur une éprouvette percée de deux trous de diamètre 1 mm. Les essais
sont réalisés sur une éprouvette en acier inox, à 850°C. L’hétérogénéité du champ de
déformation résultant de la présence des deux trous est particulièreemnt intéressante.
L’influence de la température sur les mécanismes de localisation de la déformation, ou encore
l’influence de la distance et la forme des porosités sur les mécanismes d’endommagement à
chaud sont autant de thèmes scientifiques très prometteurs pour l’avenir.
Enfin l'optimisation du cycle de vie complet d'une pièce est un grand défi. La démarche
présentée dans ce manuscrit propose l'optimisation des technologies d'assemblage par rapport
à la tenue statique en service. L'une des perspectives de ces travaux est d'optimiser le point
d'assemblage par rapport à la tenue en fatigue. À termes, ces deux objectifs pourront être
optimisés par des méthodes multi-objectifs.
208
Références
Références
[Bai08] : Y. Bai, T. Wierzbicki, A new model of metal plasticity and fracture with
pressure and Lode dependence, International Journal of Plasticity 24, pp.1071-1096, 2008.
[Bee03] : W.C.M. van Beers, J.P.C. Kleijnen, Kriging for interpolation in random
simulation, Journal of the Operational Research Society 54, pp.255–262, 2003.
209
Références
[Dix94] : L.C.W. Dixon, D.J. Mills, Effect of Rounding Errors on the Variable Metric
Method, Journals of optimization theory and application 80, pp.175-179, 1994.
[Do06] : T.T. Do, Optimisation de forme en forgeage 3D, thèse de doctorat, Mines
ParisTech, 2006.
210
Références
[Fis25] : R. Fisher, Statistical Methods for Research Worker, Olivier and Boyd, 1925.
[Gre06] : M. Grediac, F. Pierron, Applying the virtual fields method to the identification
of elasto-plastic constitutive parameters, International Journal Of Plasticity 22, pp.602-627,
211
Références
2006.
[Ham00] : V. Hamel, J.M. Roelandt, J.N. Gacel, F. Schmit, Finite element modeling of
clinch forming with automatic remeshing, Computers & Structures 77, pp.185-200, 2000.
[Han06] : A. Hana, Study and evaluation of fretting critical slip conditions by applying the
design of experiments method, Wear 261, pp.1080-1086, 2006.
[He10] : X. He, Recent development in finite element analysis of clinched joints, The
International Journal of Advanced Manufacturing Technology 48, pp.607-612, 2010.
[Hol75] : H.J. Holland, Adaptation in natural and artificial system, Ann Harbor, The
University of Michigan Press, 1975.
[Hol62] : H.J. Holland, Outline for a logical theory of adaptative systems, journal of the
association of computing machinery 9, pp.297-314, 1962.
212
Références
[Lem01] : J. Lemaitre, J.L. Chaboche, Mécanique des matériaux solides, Dunod, 2001.
213
Références
[Pil05] : M. Pillet, Plan d'expériences par la méthode Taguchi, Les Edition d'Organisation,
2005.
[Tag87] : G. Taguchi, Orthogonal arrays and linear graph, American Supplier Institute
Press, 1987.
214
Références
[Wol97] : D.H. Wolpert, W.G. Macready, No free lunch theorems for optimization, IEEE
Transactions On Evolutionary Computation 1, pp.67-82, 1997.
215
Annexe
217
Annexe
d’optimisation différentes ont été testées avec MOOPI. Dans un premier temps, il s’agissait
d’optimiser la densité de puissance dans une pièce soumise à un préchauffage avant mise en
forme. Le deuxième test consistait à atteindre une température cible en un point de la pièce.
1.6 Résultats :
• Cas 1
Pour la répartition de la densité de puissance, nous avons obtenu une très bonne
corrélation après 216 itérations (Figure A-1) :
3,00E+07
2,50E+07
Optimized power
density
2,00E+07
Power density ( W/m3 )
Computed power
density
1,50E+07
1,00E+07
5,00E+06
0,00E+00
0 0,01 0,02 0,03 0,04 0,05 0,06
Distance to the surface of the piece (m)
218
Annexe
• Cas 2
Pour l’évolution en un point de la température, nous obtenons le résultat présenté Figure
A-2 après 67 itérations.
1000
900
Profil Cible
850
800
0 20 40 60 80 100 120 140 160 180 200
Temps en secondes
Figure A-2 : Contrôle de la montée en température d'un point de la pièce - résultat de l'optimisation
La moins bonne corrélation que nous observons est principalement due au modèle direct
utilisé qui ne prend pas encore en compte la dépendance thermique de certains paramètres
magnétiques et électriques.
219
Annexe
220
Annexe
∑( ) ∑ (Tɺ ( p ) − Tɺ )
n n
Ti num ( p ) − Ti exp
2 2
num exp
i i
f ( p) = i =1
+ i =1
∑( ) ∑ (Tɺ )
n 2 n 2
Ti exp i
exp
i =1 i =1
num exp
Où p est le jeu de paramètres, T et T les observables numériques et expérimentales en
température et Tɺ représente leur dérivée temporelle.
Pour construire la loi thermo-dépendante point à point, l'hypothèse est faite que
l'absorptivité est croissante suivant la température. Il faut donc imposer une contrainte de
croissance entre les paramètres d'optimisation, une version fonction spécifique a donc été
ajoutée à MOOPI pour cette étude.
2.6 Résultats :
Les champs de température à la fin de l'impact laser et après refroidissement sont
présentés Figure A-4. Les résultats en terme d'observable sont présentés Figure A-5 et Figure
A-6. On observe une bonne corrélation à la fois pendant la phase de chauffe et pendant la
phase de refroidissement.
221
Annexe
2000
Observable numérique
1000 identifiée
500
0
0 2 4 6 8 Temps (s) 10
Figure A-5 : Profil de température associé au thermocouple situé dans la zone d'impact du
laser
250
200
Température en °C
150
100
Observable pseudo-experimentale
50
Observable numérique identifiée
0
Temps (s)
0 2 4 6 8 10
Figure A-6 : Profil de température associé au thermocouple déporté de 3mm par rapport à la
zone d'impact du laser
222
Annexe
∑ (y )
nbpo int s
2
exp
i − yinum
φPhase = i =1
∑( y )
2
exp
i
i
où y correspond aux valeurs des observables respectivement expérimentales et
numériques.
L'exploitation de cet essai complexe à plusieurs phases conduit donc à la construction d'un
problème de minimisation multi-objectifs :
Nb pahse
φ= ∑
J =1
ω J φJ
où ω J est la pondération de chaque phase.
223
Annexe
3.6 Résultats :
CONFIDENTIEL
224
Assemblage Mécanique : Stratégies d’optimisation des procédés et d’identification
des comportements mécaniques des matériaux
RESUME : Ce travail de thèse porte sur la mise en place d'une méthodologie d'optimisation des
procédés d'assemblage, ainsi que sur l'identification de paramètres de lois de comportement et
d’endommagement pour la modélisation des assemblages par déformation plastique. La première
partie du travail est axée sur la création d'une plateforme d'optimisation et d’analyse inverse
permettant d'exploiter des calculs par éléments finis. Dans le cadre du développement de cette
plateforme, un algorithme d'optimisation basé sur l'exploitation d'un méta-modèle par krigeage a été
développé. Afin de répondre aux problématiques de temps de calcul, une version parallèle de cet
algorithme exploitant les propriétés du krigeage est proposée. Cette version parallèle montre des
accélérations significatives par rapport à la version séquentielle, accélérations qui permettent un gain
de temps substantiel pour résoudre le problème d'optimisation. Le second axe de ce travail propose
une procédure d'optimisation de la tenue mécanique des assemblages par déformation plastique.
Cette optimisation est basée sur une chaîne de modélisations incluant la phase d'assemblage et de
tenue mécanique à l’arrachement. L'application de cette procédure à la technologie de clinchage
permet un gain de 13% de la tenue mécanique à l'arrachement par rapport à la solution de référence.
La troisième partie du travail est axée sur l'identification par analyse inverse de paramètres de lois de
comportement élasto-plastique et d’endommagement. L'identification, basée sur des essais de
traction, se focalise sur une loi de type élasto-plastique, couplée au modèle d'endommagement ductile
de Lemaitre. Une attention particulière est portée au type d'observables utilisées pour réaliser
l'identification. Trois types d'observables sont étudiées : la courbe force/déplacement, la mesure de
striction, et le champ de déplacement issu de mesures de champs. L'enrichissement de la base
d'observables permet de mieux définir le problème inverse en diminuant notamment les corrélations
entre paramètres et rend donc l'identification plus efficace.
ABSTRACT: This study deals with the development of a numerical strategy dedicated to
mechanical joining processes optimisation and to parameters identification by inverse analysis. The
first part is dedicated to the definition and development of an optimisation and inverse analysis
platform. This platform is adapted to finite element computations. An optimisation algorithm with a
kriging meta-model has been developed and included in the platform. In order to deal with time
consuming computation a parallel version of this algorithm, based on kriging properties, has been
developed. Significant acceleration of the parallel algorithm has been observed, leading to a
decrease in time of the resolution of the optimisation issue. In the second part of this work an
optimisation methodology has been carried out for mechanical joining processes. This procedure
enables to optimize a global simulation chain, including both the joining process and the mechanical
strength analysis of the joined component. This procedure is applied to the clinching joining process
and gives rise to an 13% increase of the mechanical strength of the component. The third part of
this work deals with a parameter identification study of an elastic-plastic law coupled with a ductile
damage model. The identification procedure is based on a tensile test. A power hardening law and a
Lemaitre model are chosen respectively for the elastic-plastic behaviour and for ductile damage.
Three different observables have been taken into account: the load/displacement curve, the necking
measurement, and displacement fields. Displacement fields are measured by full field measurement
methods. It is shown how the enrichment of the observables database improves the definition of the
inverse problem and decreases correlation issues between parameters.