AIX-MARSEILLE UNIVERSITÉ

ÉCOLE DOCTORALE no 353

SCIENCES POUR L’INGÉNIEUR : MÉCANIQUE, PHYSIQUE,
MICRO ET NANOÉLECTRONIQUE

THÈSE
pour obtenir le titre de

Docteur en Mécanique des solides

d’AIX-MARSEILLE UNIVERSITÉ
Présentée par

Hao LIU

Stratégie de raffinement automatique

de maillage et méthodes multi-grilles
locales pour le contact : application à
l’interaction mécanique Pastille-Gaine
Thèse dirigée par Frédéric LEBON
Encadrée par Isabelle RAMIERE

JURY

Mikaël BARBOTEU Professeur à l’Université de Perpignan Rapporteur

Serge DUMONT Professeur à l’Université de Nîmes Rapporteur
David DUREISSEIX Professeur à l’INSA Lyon Président du jury
Frédéric LEBON Professeur à l’Université d’Aix-Marseille Directeur de thèse
Bruno MICHEL Ingénieur-Chercheur au CEA Cadarache Examinateur
Isabelle RAMIERE Ingénieur-Chercheur au CEA Cadarache Encadrante
Remerciements

Trois ans de travail de thèse sont très enrichissants. Après toutes les périodes de doute et
d’inspiration, vient ce moment de reconnaissance.
Je voudrais tout d’abord remercier mon encadrante Isabelle Ramière et mon directeur
de thèse Frédéric Lebon pour leur accompagnement et leur soutien tout au long de cette
thèse. Ils m’ont guidé et permis de découvrir la magie de la mécanique numérique. Leur esprit
scientifique m’a beaucoup appris. C’est un honneur pour moi d’avoir travailler avec eux durant
ces trois dernières années.
Je remercie Mikaël Barboreu et Serge Dumont d’avoir accepté de relire et de rapporter
mon travail. Leurs remarques et commentaires m’ont permis d’améliorer ce mémoire. Merci à
David Dureissex de m’avoir fait l’honneur d’accepter d’être président du jury. Merci également
à Bruno Michel pour les discussions plus mécaniques ainsi que pour son aide et son intérêt
sur mon travail.
Je tiens à remercier ma hiérarchie : Renaud Masson, Bruno Collard, Carole Valot, Mireille
Bauer pour la confiance qu’ils ont mise en moi.
Je remercie mes sympathiques collègues du LSC pour m’avoir accueilli dans cet excellent
laboratoire.
Je n’oublierai jamais les super secrétaires Laure Imbert et Régine Bousquet qui se sont
toujours occupées de moi.
Merci à Jérôme Sercombe, Vincent Marelle, Thomas Helfer et Etienne Castelier pour leur
aide.
Je voudrais remercier également mes collèges LMA de m’avoir bien accueilli.
Merci à mes amis Marcelle, Jules, Coralie, Maverick, Karl, Sarah, Lei, David, Bachir, ... pour
les bons moments partagés. J’ai beaucoup aimé chaque minute que j’ai passé avec eux.
Enfin, un grand merci à ma famille pour leur amour inconditionnel et leur soutient fonda-
mental.
Stratégie de raffinement automatique de maillage et méthodes multi-grilles locales pour
le contact : application à l’interaction mécanique Pastille-Gaine

Résumé : Ce travail de thèse s’inscrit dans le cadre de l’étude de l’Interaction mécanique Pastille-Gaine
(IPG) se produisant dans les crayons combustibles des réacteurs à eau pressurisée. Ce mémoire porte
sur le développement de méthodes de raffinement de maillage permettant de simuler plus précisément
le phénomène d’IPG tout en conservant des temps de calcul et un espace mémoire acceptables pour des
études industrielles. Une stratégie de raffinement automatique basée sur la combinaison de la méthode
multi-grilles Local Defect Correction (LDC) et l’estimateur d’erreur a posteriori de type Zienkiewicz et Zhu
est proposée. Cette stratégie s’appuie sur l’erreur fournie par l’estimateur pour détecter les zones à raffiner
constituant alors les sous-grilles locales de la méthode LDC. Plusieurs critères d’arrêt sont étudiés afin de
permettre de stopper le raffinement quand la solution est suffisamment précise ou lorsque le raffinement
n’apporte plus d’amélioration à la solution globale.
Les résultats numériques obtenus sur des cas tests 2D élastiques avec discontinuité de chargement per-
mettent d’apprécier l’efficacité de la stratégie proposée.
Le raffinement automatique de maillage dans le cas de problèmes de contact unilatéral est ensuite abordé.
La stratégie proposée dans ce travail s’étend aisément au raffinement multi-corps à condition d’appliquer
l’estimateur d’erreur sur chacun des corps séparément. Un post-traitement est cependant souvent néces-
saire pour garantir la conformité des zones de raffinement vis-à-vis des frontières de contact. Une variété
de tests numériques de contact (avec ou sans frottement, avec ou sans jeu initial) entre solides élastiques
confirme l’efficacité et la généricité de la stratégie proposée.

Mots clés : Raffinement adaptatif de maillage, Méthodes multi-grilles locales, Estimateur d’erreur a
posteriori, Interaction mécanique Pastille-Gaine, Contact unilatéral, Loi de frottement de Coulomb

Automatic mesh refinement and local multigrid methods for contact problems :
application to the Pellet-Cladding mechanical Interaction

Abstract : This Ph.D. work takes place within the framework of studies on Pellet-Cladding mechanical
Interaction (PCI) which occurs in the fuel rods of pressurized water reactor. This manuscript focuses on
automatic mesh refinement to simulate more accurately this phenomena while maintaining acceptable
computational time and memory space for industrial calculations. An automatic mesh refinement strategy
based on the combination of the Local Defect Correction multigrid method (LDC) with the Zienkiewicz
and Zhu a posteriori error estimator is proposed. The estimated error is used to detect the zones to be
refined, where the local subgrids of the LDC method are generated. Several stopping criteria are studied
to end the refinement process when the solution is accurate enough or when the refinement does not
improve the global solution accuracy anymore.
Numerical results for elastic 2D test cases with pressure discontinuity shows the efficiency of the proposed
strategy.
The automatic mesh refinement in case of unilateral contact problems is then considered. The strategy
previously introduced can be easily adapted to the multibody refinement by estimating solution error
on each body separately. Post-processing is often necessary to ensure the conformity of the refined areas
regarding the contact boundaries. A variety of numerical experiments with elastic contact (with or without
friction, with or without an initial gap) confirms the efficiency and adaptability of the proposed strategy.

Keywords : Automatic mesh refinement, Local multigrid method, a posteriori error estimator, Pellet-
Cladding mechanical Interaction, Unilateral contact, Coulomb’s friction law
Table des matières

Introduction 9

1 Contexte industriel 11
1.1 Réacteur à Eau Pressurisée (REP) . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Crayon combustible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1 Pastille de combustible . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.2 Gaine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3 Comportement du crayon combustible . . . . . . . . . . . . . . . . . . . . . . 16
1.3.1 Dilatation thermique . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.2 Densification et gonflement du combustible . . . . . . . . . . . . . . . 18
1.3.3 Fluage de la gaine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Interaction mécanique pastille-gaine . . . . . . . . . . . . . . . . . . . . . . . 18
1.5 Objectif du travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2 Raffinement de maillage en mécanique des solides 23

2.1 Préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Problème élastostatique linéaire . . . . . . . . . . . . . . . . . . . . . . 25
2.1.2 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.3 Approximation par la méthode des éléments finis . . . . . . . . . . . . 28
2.1.4 Erreur de discrétisation . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Raffinement de maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.1 Méthodes adaptatives . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.2 Méthodes multi-grilles locales . . . . . . . . . . . . . . . . . . . . . . . 35
2.3 Estimateurs d’erreur a posteriori . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.2 Critères de comparaison des estimateurs . . . . . . . . . . . . . . . . . 40
2.3.3 Estimateur basé sur le résidu d’équilibre . . . . . . . . . . . . . . . . . 40
2.3.4 Estimateur basé sur la loi de comportement . . . . . . . . . . . . . . . 45
2.3.5 Estimateur basé sur le lissage des contraintes . . . . . . . . . . . . . . 47

3 Stratégie de raffinement automatique de maillage 51

3.1 Méthode “Local defect correction” . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.1.1 Formulation de la méthode LDC . . . . . . . . . . . . . . . . . . . . . . 54
8 Table des matières

3.1.2 Application la méthode LDC dans le cadre de la méthode des éléments

finis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2 Combinaison de la méthode LDC et de l’estimateur ZZ . . . . . . . . . . . . . 60
3.3 Présentation des cas tests 2D issus de la simulation de l’IPG . . . . . . . . . . 63
3.3.1 Mise en diabolo : cas 2D(r,z) . . . . . . . . . . . . . . . . . . . . . . . 63
3.3.2 Fissuration de la pastille : cas 2D(r,θ) . . . . . . . . . . . . . . . . . . . 64
3.4 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.1 Choix numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.2 Erreur relative en norme énergie . . . . . . . . . . . . . . . . . . . . . 65
3.4.3 Erreur absolue en norme infinie . . . . . . . . . . . . . . . . . . . . . . 77

4 Application de la stratégie de raffinement automatique au problème de contact 83

4.1 Problème de contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.1.1 Problème de contact unilatéral . . . . . . . . . . . . . . . . . . . . . . 86
4.1.2 Problème de contact avec frottement . . . . . . . . . . . . . . . . . . . 88
4.2 Quelques méthodes pour résoudre le problème de contact avec ou sans frottement 91
4.2.1 Méthode de pénalisation . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.2.2 Méthode des multiplicateurs de Lagrange . . . . . . . . . . . . . . . . 92
4.2.3 Méthode du Lagrangien augmenté . . . . . . . . . . . . . . . . . . . . 92
4.2.4 Méthode de point fixe . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.2.5 Algorithmes de résolutions . . . . . . . . . . . . . . . . . . . . . . . . 93
4.3 Stratégie de raffinement automatique du maillage adaptée au problème de
contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.4 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.4.1 Présentation du cas test . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.4.2 Illustration de la détection des zones à raffiner . . . . . . . . . . . . . 97
4.4.3 Résultats numérique pour un cas test sans jeu ni frottement . . . . . . 98
4.4.4 Résultats numériques pour un cas test avec jeu sans frottement . . . . 105
4.4.5 Résultats numériques pour un cas test avec jeu avec frottement . . . . 106
4.4.6 Contact avec frontière courbe . . . . . . . . . . . . . . . . . . . . . . . 110

Conclusion et Perspectives 117

Bibliographie 119
Introduction

De la recherche à l’industrie, la simulation numérique par la méthode des éléments finis est
devenue un outil essentiel permettant de modéliser des phénomènes physiques divers. Cette
méthode associée au développement de méthodes numériques performantes de résolution
permet de simuler et de comprendre des phénomènes et le comportement d’objets de plus
en plus complexes. La précision de la solution par éléments finis est directement liée à la
taille des éléments du maillage. La conséquence est que les modèles numériques présentent
potentiellement un très grand nombre de degrés de liberté.
En vue d’atteindre des tailles de maille localement fines sans avoir à raffiner uniformément
le maillage pour des raisons évidentes de temps de calcul et d’espace mémoire, différentes
méthodes de raffinement de maillage ont été développées, en particulier les méthodes de
raffinement adaptatif et les méthodes multi-grilles locales. Les méthodes de raffinement adap-
tatif ont pour but de modifier le maillage initial dans des zones d’intérêt telles que singularité
ou concentration de contraintes afin d’obtenir un maillage optimal en terme de taille ou
degré d’élément. Les méthodes multi-grilles locales consistent quant à elles à ajouter des
sous-niveaux avec un pas de maillage de plus en plus fin dans ces zones d’intérêt et de re-
lier les solutions de chaque niveau par un processus itératif. Afin de permettre de détecter
automatiquement les zones d’intérêt, les méthodes de raffinement de maillage sont souvent
couplées avec un estimateur d’erreur a posteriori.
En mécanique des solides, le phénomène de contact entre solides déformables est un pro-
blème largement répandu et étudié. La résolution précise des tels problèmes par la méthode
des éléments finis requiert des mailles suffisamment fines autour des zones de contact et
est donc souvent assez coûteuse. À titre d’exemple, dans le cadre de l’étude de l’Interaction
mécanique Pastille combustible-Gaine (IPG) dans les Réacteurs à Eau Pressurisée (REP), la
représentation précise des phénomènes mis en jeu nécessite des mailles fines de l’ordre du
micron autour de la zone de contact pastille-gaine pour une structure de l’ordre du centi-
mètre. Au-delà de la formulation du problème de contact et de sa résolution mathématique,
beaucoup de travaux de recherche et développement s’intéressent désormais à l’optimisation
locale du maillage de calcul.
Cette thèse porte sur le développement de stratégies de raffinement de maillage permettant
de simuler plus précisément le phénomène d’IPG avec des temps de calcul et un espace
mémoire acceptables. En effet, la gaine est la première barrière de confinement des produits
de fission. Son intégrité et son étanchéité doivent être garanties, notamment vis-à-vis du
phénomène de rupture. Le phénomène d’IPG doit donc être bien compris et simulé dans la
10 INTRODUCTION

démarche de sûreté de l’élément combustible des REP.

Cette thèse a été réalisée dans le cadre du projet SICOM (SImulation COMbustible) dédié
au développement de la plate-forme logicielle PLEIADES pour la simulation du comportement
des éléments combustibles sous irradiation. Ces travaux sont cofinancés par le CEA (Com-
missariat à l’Énergie Atomique et aux Énergies Alternatives), EDF (Électricité de France) et
AREVA.

Plan de ce mémoire
Ce mémoire de thèse se compose de quatre chapitres.
Le chapitre 1 détaille le contexte industriel. Nous y présentons les différentes phénomènes
mécaniques se produisant à l’intérieur du combustible d’un réacteur à eau pressurisée.
Dans le chapitre 2, les différentes techniques de raffinement de maillage pour un problème
élastostatique linéaire et les méthodes d’estimation d’erreur a posteriori permettant de détecter
automatiquement les zones à raffiner sont présentées.
Nous proposons dans le chapitre 3 une stratégie de raffinement automatique de maillage
basée sur le couplage de la méthode Local Defect Correction et l’estimateur d’erreur a posteriori
de Zienkiewicz et Zhu. Cette stratégie est vérifiée sur des cas test 2D élastiques représentatifs
de l’IPG.
Enfin, dans le chapitre 4, nous étendons la stratégie proposée à des problèmes de contact
avec ou sans frottement. Nous vérifions notamment qu’un raffinement multi-corps s’effectue
bien autour des zones de contact.
Notre objectif est de proposer une stratégie de raffinement local automatique unique et gé-
nérique permettant de raffiner différents types de problèmes comportant des singularités
locales, notamment des problèmes de contact, ce qui serait novateur vis-à-vis de la littérature
existante.
Chapitre 1

Contexte industriel

Dans ce chapitre, nous présentons le contexte industriel de la thèse. Il concerne les phéno-
mènes mécaniques se produisant à l’intérieur du combustible d’un réacteur à eau pressurisée
et qui conduisent, sous l’effet de la température et de l’irradiation, à l’interaction mécanique
entre la pastille et la gaine de combustible.

Sommaire
1.1 Réacteur à Eau Pressurisée (REP) . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Crayon combustible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1 Pastille de combustible . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.2 Gaine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3 Comportement du crayon combustible . . . . . . . . . . . . . . . . . . . . . . 16
1.3.1 Dilatation thermique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.2 Densification et gonflement du combustible . . . . . . . . . . . . . . . . 18
1.3.3 Fluage de la gaine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Interaction mécanique pastille-gaine . . . . . . . . . . . . . . . . . . . . . . . 18
1.5 Objectif du travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.1. Réacteur à Eau Pressurisée (REP) 13

Les réacteurs à eau légère et en particulier les réacteurs à eau pressurisée (REP) dominent
le marché éléctronucléaire international depuis plusieurs décennies. En effet sur les 430
réacteurs en opération à travers le monde, 2/3 sont des REP, et sur les 70 réacteurs en
constructions, 3/4 sont également des REP.
Avec ses 58 réacteurs nucléaires répartis dans 19 centrales nucléaires produisant près
de 80% de l’électricité nationale, la France est un pays fortement axé autour de l’énergie
nucléaire. Tous les réacteurs nucléaires en service en France sont des REP de puissance
900 MWe à 1450 MWe (MégaWatt électrique, unité de puissance électrique). Nous rappelons
ici brièvement le fonctionnement de ce type de réacteur.

1.1 Réacteur à Eau Pressurisée (REP)

Un réacteur nucléaire à eau pressurisée (type REP) est composé de 3 circuits (voir fi-
gure 1.1) :

F I G U R E 1.1 – Centrale nucléaire à eau pressurisée

• Le circuit primaire
Dans le coeur du réacteur, la fission nucléaire produit de la chaleur. De l’eau, jouant le
rôle de modérateur et de caloporteur, circule dans le circuit primaire à la pression de
14 Chapitre 1. Contexte industriel

15,5 MPa. Cette eau, à la température d’environ 290 ◦ C à l’entrée du coeur, est chauffée
au contact des crayons combustibles à environ 324 ◦ C en sortie du coeur.
• Le circuit secondaire
Dans le générateur de vapeur, l’eau du circuit secondaire est chauffée par le circuit
primaire au contact des tubes et se transforme en vapeur. Cette vapeur fait tourner une
turbine couplée à un alternateur qui produit de l’électricité.
• Le circuit de refroidissement
L’eau de ce troisième circuit refroidit le circuit secondaire à travers un échangeur et
condenseur.
Vue la radioactivité des produits de fission nucléaire, les exigences de sûreté sont fortes pour
éviter l’éventuelle dispersion d’éléments radioactifs à l’extérieur de la centrale. Il existe ainsi
3 barrières de confinement :
• Première barrière : un gainage métallique qui enferme le combustible
• Seconde barrière : une enveloppe métallique autour du circuit primaire (tube de Géné-
rateur de vapeur, cuve de réacteur, ...)
• Troisième barrière : une enceinte de confinement en béton de forte épaisseur

1.2 Crayon combustible

Le crayon combustible est un élément de base du réacteur REP. C’est l’endroit où la cha-
leur est produite. Un crayon combustible de REP est constitué essentiellement d’une gaine
contenant environ 300 pastilles de combustible maintenues par un ressort (voir figure 1.2). Le
crayon, pressurisé à l’hélium entre 2 MPa et 2,5 MPa à la fabrication, est fermé par deux bou-
chons soudés à ses extrémités. L’hélium est choisi car il a la meilleure conductivité thermique
parmi les gaz nobles. Le jeu radial initial entre les pastilles et la gaine est d’environ 85 µm
en vue de faciliter l’introduction des pastilles dans la gaine. Nous détaillons dans la suite les
principales caractéristiques des pastilles de combustible et de la gaine ainsi que leurs condi-
tions de fonctionnement. Pour plus de détails, le lecteur peut se reporter à [Bailly et al. 1997]
et [Guerin et al. 2013].
Bouchon Bouchon
supérieur Ressort de maintien Pastilles inférieur

17,45mm 3660mm 17,45mm

3852mm

F I G U R E 1.2 – Dimension d’un crayon combustible pour un REP de puissance de 900 MWe
1.2. Crayon combustible 15

1.2.1 Pastille de combustible

Le combustible des REP se présente sous forme de pastilles cylindriques mesurant 8,2 mm
de diamètre pour 13,5 mm de hauteur (voir figure 1.3). Le matériau est du dioxyde d’uranium
(UO2 ) ou de l’oxyde mixte d’uranium et de plutonium (MOX). La pastille est une céramique
fragile avec une contrainte de rupture comprise entre 100 MPa et 150 MPa. Après fabrication,
la pastille a une porosité comprise entre 5 et 6%.
Oz axe du crayon combustible
6 mm
Chanfrein Plan inter-pastille

13, 5 mm
Plan median-pastille

Evidement

Plan inter-pastille
8, 2 mm

(a) Pastilles de combustible (b) Dimensions caractéristiques des pastilles de

combustible

F I G U R E 1.3 – Pastilles de combustible utilisées dans les REP

La pastille de combustible présente un évidement hémisphérique et un chanfrein à chacune

de ses extrémités. L’évidement a pour but de compenser l’excès de dilatation du centre de la
pastille et le chanfrein facilite son introduction dans la gaine.

1.2.2 Gaine
La gaine est un tube de diamètre extérieur valant 9,5 mm et d’épaisseur 0.57 mm. Il s’agit
de la première barrière de confinement des produits de fission.
Les matériaux utilisés pour la gaine des REP actuellement en fonctionnement en France et
à l’international sont des alliages à base de zirconium (Zircaloy-4,M5
R
, ZIRLO TM , Optimized
ZIRLOTM ). Leurs propriétés répondent à certains besoins :
• une très faible absorption des neutrons thermiques (utiles à la réaction en chaine),
• une bonne conductivité thermique en vue d’évacuer vers l’eau la chaleur produite dans
le crayon,
• une bonne résistance mécanique,
• une bonne résistance à la corrosion vis-à-vis de l’eau à haute température et du com-
bustible.
16 Chapitre 1. Contexte industriel

1.3 Comportement du crayon combustible

Le comportement du crayon combustible dans un réacteur nucléaire est principalement
piloté par la haute température et l’irradiation. La connaissance des phénomènes mécaniques
concernant la pastille et la gaine ainsi que leurs interactions est très importante pour s’assurer
du bon fonctionnement du crayon et de l’assemblage combustible.
Le schéma 1.4 décrit les principaux phénomènes se produisant dans un crayon REP du
début de sa vie jusqu’au 3e cycle 1 de fonctionnement en réacteur. Dans cette section, nous
allons uniquement présenter les phénomènes relatifs à l’interaction mécanique pastille-gaine.

A) Combustible neuf (hors réacteur) B) Montée en puissance

- 3 MPa à froid dans le jeu combustible/gaine - Dialatation différentielle
- Pression externe 15,5MPa - Fracturation de la pastille
- Diminution du jeu

D) 3e cycle C) 1er et 2 e cycle

- Corrosion et hydruration de la gaine - Fluage de la gaine vers le combustible
- Réarrangement des fragments - Densification
- La gaine exerce une contrainte sur le combustible - Fermeture du jeu
- Apparition d'une zone poreuse - Effet diabolo

F I G U R E 1.4 – Comportement générale d’un crayon REP [Bailly et al. 1997]

1.3.1 Dilatation thermique

Dans les REP, la température dans le coeur d’une pastille peut atteindre plus de 1000 ◦ C
en fonctionnement normal. Le premier phénomène présenté ici est la dilatation thermique de
la pastille et de la gaine. Du fait que les coefficients de dilatation et les températures de la
pastille et de la gaine sont différentes, leurs dilatations sont également différentes. Bien que le
coefficient de dilatation de la gaine soit plus important que celui de la pastille, la température
dans la pastille est beaucoup plus élevée que celle dans la gaine (voir figure 1.5). La dilatation
de la pastille est finalement plus importante que celle de la gaine. Par conséquent, le jeu radial
initial d’environ 85 µm diminue jusqu’à environ 60 µm en raison des dilatations thermiques
différentielles de la pastille et de la gaine.
1. un cycle correspond au temps entre deux arrêts du réacteur pour rechargement du coeur et maintenance,
soit environ 12 à 18 mois suivant les modes de gestion. Un cycle comprend donc une montée de puissance
(démarrage du réacteur), un maintien de puissance plus ou moins constant (fonctionnement du réacteur) et un
retour à froid (arrêt du réacteur)
1.3. Comportement du crayon combustible 17

F I G U R E 1.5 – Exemple du profil radial de la température dans le crayon combustible

[Guerin et al. 2013]

La conductivité thermique de la pastille étant faible et la puissance importante, un fort

gradient radial de température est présent dans celle-ci. La figure 1.5 présente un exemple du
profil radial de température dans le crayon combustible. La différence de dilatation thermique
à l’intérieur de la pastille sous l’effet du fort gradient de température radial provoque de fortes
contraintes dans la pastille. Dès la fin de la première montée en puissance, les pastilles sont
systématiquement fissurées (voir figure 1.6). L’orientation des fissures est majoritairement
radiale.
Sous l’effet des contraintes thermiques, les déformations dans le plan inter-pastilles sont
plus importantes que celles dans le plan médian-pastille. La pastille se déforme sous la forme
de diabolo (voir figure 1.7).

F I G U R E 1.6 – Macrographie d’une pastille REP fis- F I G U R E 1.7 – Mise en diabolo de la

surée par le gradient thermique pastille
18 Chapitre 1. Contexte industriel

1.3.2 Densification et gonflement du combustible

En début de vie, l’irradiation conduit au phénomène de densification. La baisse de la
porosité diminue le volume de la pastille et retarde légèrement l’instant du premier contact
pastille-gaine.
Au cours de l’irradiation, la production de produits de fission au sein du combustible
provoque une déformation macroscopique, appelée gonflement de la pastille. Ce gonflement
dépend donc de la quantité et de l’état physico-chimique des produits de fission.

1.3.3 Fluage de la gaine

Le fluage est le phénomène physique qui provoque progressivement une déformation ir-
réversible sous charge constante inférieure à la limite d’élasticité du matériau. Il existe deux
types de fluage pour la gaine : le fluage d’irradiation et le fluage thermique [Bailly et al. 1997].
Le fluage d’irradiation décrit la déformation sous l’effet du flux neutronique. Le flux neutro-
nique participe à la mobilité des défauts et accroît la vitesse de fluage. La vitesse du fluage
thermique dépend de la température.
En fonctionnement normal, à faible contrainte (inférieure à 100 MPa), le fluage d’irra-
diation contribue à près de 80% de la déformation totale par fluage. À forte contrainte, il
représente environ 20% de la déformation totale par fluage.
La pression interne dans le crayon est initialement de 2,5 MPa et augmente jusqu’à environ
6 MPa en début d’irradiation. Au cours de la vie, dû au relâchement des gaz de fission, cette
pression continue d’augmenter mais ne dépasse pas celle du fluide caloporteur. La gaine subit
donc une contrainte circonférentielle de compression. Sous l’effet de l’irradiation, le diamètre
de la gaine a donc tendance à diminuer.

1.4 Interaction mécanique pastille-gaine

Au cours du premier cycle, le jeu initial entre la pastille et la gaine diminue du fait des
phénomènes précedément cités : gonflement et mise en diabolo de pastille, et fluage de la
gaine. Au début du deuxième cycle d’irradiation, le jeu commence à se fermer sur certaines
zones du crayon et un contact entre la pastille et la gaine se présente au niveau des plans inter-
pastilles. Ce phénomène est appelé l’Interaction Pastille-Gaine (IPG). De plus, la fissuration de
la pastille introduit un contact discontinu entre la pastille et la gaine. Ces deux phénomènes
induisent une concentration locale importante de contraintes au niveau du contact pastille-
gaine, au droit des fissures.
La fermeture du jeu entre la pastille et la gaine modifie les contraintes auxquelles la
gaine est soumise. Avant la fermeture du jeu, la contrainte principale dans la gaine est une
contrainte de compression due au fluide calopoteur. Lorsque la fermeture du jeu a eu lieu, la
valeur absolue de la contrainte de compression décroît jusqu’à devenir nulle et devient ensuite
une contrainte de traction à cause des contraintes exercées par la pastille. Après plusieurs
cycles d’irradiation, la gaine présente une forme de bambou due à l’influence de la mise
en diabolo des pastilles (voir figure 1.8). Au niveau du plan inter-pastilles, un déplacement
1.4. Interaction mécanique pastille-gaine 19

plus important est observé. La figure 1.9 présente un exemple du profil axial du diamètre du
crayon après 2 cycles annuels et montre que ce déplacement reste modéré et est inférieur 20
µm en fonctionnement normal.

mesures expérimentales après irradiation

niveaux inter-pastilles

Diamètre extérieur gaine (µm)

Profil axial du diamètre du crayon(mm)

F I G U R E 1.8 – Représenta- F I G U R E 1.9 – Profil axial du diamètre du crayon après 2

tion schématique de la dé- cycles annuels [Guerin et al. 2013]
formation en bambou de la
gaine

Dans les REP, le phénomène d’IPG ne provoque pas la rupture de la gaine en fonction-
nement normal. Cependant, comme la gaine est la première barrière de confinement du
combustible, son intégrité et son étanchéité doivent être garanties non seulement en régime
normal mais également en régime incidentel de niveau 2 (probabilité de 10−2 par réacteur et
par an) selon l’échelle INES (de l’anglais International Nuclear Event Scale).
De nombreuses études sur ce phénomène ont été réalisées en vue de qualifier les li-
mites technologiques à rupture de la gaine et contraintes exercées sur le rayon interne de la
gaine. En particulier, un grand nombre d’essais de rampe [Garlick 1973, Alberman et al. 1997,
Mougel et al. 2004], simulant une situation incidentelle de niveau 2, a été effectué dans des
Réacteurs d’Essai des Matériaux (MTR acronyme Anglais pour Material Testing Reactor) afin
de déterminer la limite à rupture. Les essais expérimentaux sont cependant généralement
complexes et coûteux.
Un type d’étude complémentaire repose sur la modélisation des mécanismes condui-
sant à l’IPG [Jernkvist 1995, Michel et al. 2008b, Michel et al. 2008a, Marchal et al. 2009,
Sercombe et al. 2012, Michel et al. 2013]. Les objectifs de la modélisation sont également
de déterminer les limites à rupture du crayon ainsi que de prédire le comportement du crayon
en situation incidentelle.
En vue de simuler le comportement des éléments combustibles sous irradiation, la plate-
forme PLEIADES (Plateforme Logicielle pour les Elements Irradiés dans les Assemblages en
20 Chapitre 1. Contexte industriel

Démonstration, en Expérimentation ou en Service) a été développée par le CEA et ses par-

tenaires depuis 2003. Dans cette plate-forme, plusieurs applications sont développées pour
s’adapter aux différentes filières de réacteur étudiées. L’ensemble de ces applications par-
tage des composants communs, des lois et des modèles de comportement. Au sein de cette
plate-forme, l’application ALCYONE est destinée à la filière REP. Dans le but d’obtenir des
contraintes et des déformations locales maximales, les études de l’IPG utilisent principalement
le modèle 3D d’ALCYONE [Michel et al. 2008b, Michel et al. 2008a, Sercombe et al. 2012,
Michel et al. 2013, Sercombe et al. 2013]. En supposant que le combustible se fissure de façon
régulière (voir [Nonon et al. 2004]) et en raison de symétries, seul un fragment correspondant
à 1/32 d’une pastille et de la gaine qui l’entoure est modélisé (voir figure 1.10).

F I G U R E 1.10 – Modélisation d’1/32 d’une pastille et de la gaine

Les études expérimentales [Mougel et al. 2004] et numériques [Michel et al. 2008a,
Sercombe et al. 2013] ont confirmé que les contraintes et les déformations de la gaine sont
maximales au niveau du plan inter-pastilles en face des fractures radiales de la pastille (voir
figure 1.11). Ce point est appelé point triple.

1.5 Objectif du travail

Afin de pouvoir obtenir de manière précise les contraintes au niveau du point triple,
la taille de maille de calcul doit être localement de l’ordre du micron pour une structure
d’ordre du centimètre [Sercombe et al. 2013]. Cependant, les temps de calcul et ressources
mémoire associés à des maillages raffinés localement, sans être trop dégénérés, deviennent
vite prohibitifs pour les codes de R&D et à fortiori pour les codes industriels.
Une alternative, basée sur des méthodes multi-grilles locales permettant d’atteindre cette
taille de maille avec des temps de calculs et espaces mémoires acceptables, a été étudiée par
Barbié et al. [Barbié 2013, Barbié et al. 2014, Barbié et al. 2015]. Les auteurs ont étudié la
1.5. Objectif du travail 21

F I G U R E 1.11 – Concentration de contraintes au point triple [Sercombe et al. 2013]

performance d’une méthode multi-grilles locale de type Local Defect Correction (LDC) pour
des problèmes linéaires et des problèmes non-linéaires avec loi de fluage de Norton. Via
l’estimateur d’erreur a posteriori de Zienkiewicz et Zhu [Zienkiewicz & Zhu 1987] (ZZ), les
auteurs ont choisi de raffiner les éléments du maillage dont l’indicateur d’erreur était supérieur
à α% de l’erreur maximale où α désigne un paramètre. Cette stratégie de raffinement est
efficace et simple à mettre en oeuvre, mais elle ne permet pas d’arrêter la génération de
sous-niveaux. De plus, le paramètre α dépend du problème et est donc difficile à connaître a
priori.
Le but de ce travail est tout d’abord de développer une stratégie de raffinement auto-
matique, toujours basée sur la combinaison de la méthode LDC et de l’estimateur ZZ, mais
permettant de générer de manière optimale et générique les sous-grilles locales nécessaires
au calcul précis et ce, indépendamment de la connaissance du problème traité.
Le deuxième objectif de la thèse est de s’intéresser au raffinement multi-grilles local multi-
corps (pastille-gaine ici) et d’y étendre la stratégie de raffinement automatique. En particulier,
le raffinement doit pouvoir traiter des situations de contact entre les corps.
Chapitre 2

Raffinement de maillage en
mécanique des solides

Dans ce chapitre, nous commençons par exposer le problème modèle élastostatique li-
néaire sur lequel un rappel de la méthode des éléments finis est formulé. Ensuite, nous pré-
sentons les différentes techniques de raffinement de maillage pour ce type de problème ainsi
que les méthodes d’estimation d’erreur a posteriori permettant de détecter automatiquement
les zones à raffiner.

Sommaire
2.1 Préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Problème élastostatique linéaire . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.2 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.3 Approximation par la méthode des éléments finis . . . . . . . . . . . . . 28
2.1.4 Erreur de discrétisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Raffinement de maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.1 Méthodes adaptatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.2 Méthodes multi-grilles locales . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3 Estimateurs d’erreur a posteriori . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.2 Critères de comparaison des estimateurs . . . . . . . . . . . . . . . . . . 40
2.3.3 Estimateur basé sur le résidu d’équilibre . . . . . . . . . . . . . . . . . . 40
2.3.4 Estimateur basé sur la loi de comportement . . . . . . . . . . . . . . . . 45
2.3.5 Estimateur basé sur le lissage des contraintes . . . . . . . . . . . . . . . 47
2.1. Préliminaires 25

2.1 Préliminaires
Dans cette section, nous faisons un bref rappel du problème d’élasticité linéaire statique
et de la méthode des éléments finis utilisée pour le résoudre.

2.1.1 Problème élastostatique linéaire

On considère un solide occupant le domaine borné Ω ⊂ Rm , m = 2 ou 3 avec une frontière
suffisamment régulière Γ = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅.
On suppose que le matériau est élastique linéaire. Sous l’hypothèse des petites perturba-
tions et en statique, l’équation d’équilibre est une relation entre le tenseur des contraintes
d’ordre 2, σ, et la densité volumique de force f :

− div σ = f dans Ω. (2.1.1)

La relation entre le tenseur des contraintes σ et le tenseur des déformations linéarisé

d’ordre 2, ε, est exprimée par la loi de comportement du matériau (ou équation constitutive) :

σ = Cε (2.1.2)

avec C le tenseur d’élasticité d’ordre 4 qui admet les symétries suivantes :

Cijkl = Cjikl = Cklij (2.1.3)

et en notation d’Einstein :
(Cε)ij = Cijkl εkl (2.1.4)

Le tenseur des déformations ε est relié au déplacement u par :

1
gradu + gradT u


ε= (2.1.5)
2
Ce problème admet une solution unique sous les conditions aux limites suivantes :
— Conditions aux limites de Dirichlet : u = u0 sur ΓD avec u0 déplacement imposé
sur le bord ΓD , avec ΓD 6= ∅.
— Conditions aux limites de Neumann : σn = FN sur ΓN avec (σn)i = σij nj , n
vecteur unitaire normal dirigé vers l’extérieur de Ω et FN la force exercée sur ΓN .
Le problème élastostatique linéaire complet peut donc s’exprimer sous la forme :



 −div σ = f dans Ω
σ = Cε dans Ω




1

T


ε= grad u + grad u dans Ω (2.1.6)

 2



 u = u0 sur ΓD

 σn = FN sur ΓN
26 Chapitre 2. Raffinement de maillage en mécanique des solides

Champs de déplacement cinématiquement admissibles et de contraintes statiquement

admissibles :
Nous rappelons ici les définitions des champs admissibles [Bonnet & Frangi 2007] :
— un champ continu et régulier de déplacement est cinématiquement admissible (CA)
s’il satisfait les conditions aux limites de Dirichlet sur le bord ΓD :

uCA = {v; v continu et régulier sur Ω et v = u0 sur ΓD } (2.1.7)

— un champ régulier de contraintes est statiquement admissible (SA) s’il satisfait l’équa-
tion d’équilibre (2.1.1) et les conditions aux limites de Neumann sur le bord ΓN :

σSA = {τ ; −div τ = f dans Ω, τ n = FN sur ΓN } (2.1.8)

Remarque : la solution du problème (2.1.6) est un champ de déplacement cinématique-

ment admissible et un champ régulier de contrainte statiquement admissible associé.

2.1.2 Formulation variationnelle

Par souci de simplification, on choisira par la suite des conditions aux limites de Dirichlet
homogène : u0 = 0. L’extension des calculs et des démonstrations au cas non homogène ne
pose pas de problème (cf. remarque page 27) et conduirait à alourdir la présentation.
On note H 1 (Ω) l’espace de Sobolev défini par :

H 1 (Ω) = {v ∈ L2 (Ω); ∂i v ∈ L2 (Ω), i = 1, · · · , m}

où L2 est l’ensemble des fonctions de carré intégrable :

Z 1/2
2
L (Ω) = {v mesurable, v · vdω < ∞}
Ω

avec m la dimension de l’espace.

Soit V 0 l’espace vectoriel des champs de déplacements admissibles :

V 0 = {v ∈ (H 1 (Ω))m ; v(x) = 0, x ∈ ΓD }.

Soit une fonction test v ∈ V 0 . En multipliant l’équation équilibre (2.1.1) par v et en

l’intégrant sur le domaine Ω, on obtient :
Z Z
− divσ · v dΩ = f · v dΩ
Ω Ω

En supposant v, σ suffisamment régulières et en intégrant par parties, on a :

Z Z Z
σ : grad v dΩ − (σn) · v dΓ = f · v dΩ
Ω Γ Ω

avec
∂vi
σ : grad v = σij .
∂xj
2.1. Préliminaires 27

En utilisant v ∈ V 0 sur ΓD , on a :
Z Z Z
σ : grad v dΩ − (σn) · v dΓ = f · v dΩ
Ω ΓN Ω

Les propriétés de symétrie du tenseur des contraintes (σij = σji ) nous donnent :

1
Z Z
gradv + gradT v


σ : grad v dΩ = σ:
Ω 2
ZΩ Z
= σ : ε(v) dΩ = (Cε(u)) : ε(v) dΩ (2.1.9)
Ω Ω

On obtient finalement la formulation variationnelle du problème (2.1.6) en utilisant la condi-

tion aux limites sur ΓN du problème (2.1.6) :
Z Z Z
(Cε(u)) : ε(v) dΩ = f · v dΩ + FN · v dΓ (2.1.10)
Ω Ω ΓN

Pour simplifier l’expression et être plus général, on utilise les notations suivantes :
Z
a(u, v) := (Cε(u)) : ε(v) dΩ
ZΩ Z
l(v) := f · v dΩ + FN · v dΓ
Ω ΓN

avec a un opérateur bilinéaire symétrique, l un opérateur linéaire.

La formulation variationnelle (2.1.10) peut s’écrire sous la forme :

Trouver u ∈ V 0 tel que

(
(2.1.11)
a(u, v) = l(v) ∀v ∈ V 0

Le problème (2.1.6) est équivalent au problème (2.1.11).

Pour rappel, soit :
— V 0 un espace de Hilbert muni de son produit scalaire noté h·, ·i, de norme associée
notée k · k,
— a : V 0 ×V 0 → R une forme bilinéaire continue et coercive (ou V-elliptique) sur V 0 ×V 0 ,
c’est-à-dire,
— ∃M > 0, ∀(u, v) ∈ V 0 × V 0 , |a(u, v)| ≤ M kukkvk
— ∃ α > 0, ∀u ∈ V 0 , a(u, u) ≥ αkuk2
— l une forme linéaire continue c’est-à-dire ∃L > 0, ∀v ∈ V 0 , |l(v)| ≤ Lkvk
Alors, selon le théorème de Lax-Milgram [Lax & Milgram 1954] , il existe une unique
solution u ∈ V 0 au problème (2.1.11).
Remarque : on peut utiliser par exemple la technique de relèvement pour traiter le pro-
blème avec condition aux limites de Dirichlet non homogène (u0 6= 0). Si u0 est suffisament
régulier, d’après le théorème de trace, on peut trouver un relèvement ru de cette condition
aux limites tel que γ0 ru = u0 où γ0 est l’opérateur de trace. On cherche maintenant u sous la
forme u = ũ + ru où ũ est la solution de l’équation (2.1.12) avec des conditions aux limites
de Dirichlet homogène.
28 Chapitre 2. Raffinement de maillage en mécanique des solides



 −div Cε(ũ) = f + div Cε(ru ) dans Ω
 1
grad ũ + gradT ũ


ε(ũ) = dans Ω

2 (2.1.12)


 ũ = 0 sur ΓD

σn = FN sur ΓN


2.1.3 Approximation par la méthode des éléments finis

Soit Vh0 ⊂ V 0 un sous-espace de dimension finie, tel que Vh0 → V 0 lorsque h → 0. On
cherche ici une approximation uh ∈ Vh0 de la solution u ∈ V 0 . Le problème approché s’écrit :

Trouver uh ∈ Vh0 tel que

(
(2.1.13)
a(uh , vh ) = l(vh ) ∀vh ∈ Vh0

À titre illustratif, nous allons présenter la méthode des éléments finis de Lagrange en 2D
pour des éléments triangulaires. Pour d’autres types d’éléments, le lecteur pourra se reporter
à [Raviart & Thomas 1998]. Soit Ω un ouvert polygonal de R2 , on se donne une triangulation
Th de Ω, c’est-à-dire un ensemble de triangles (Ki )i=1,··· ,nk tel que Ki ⊂ Ω̄, Ω̄ = ∪ni=1k
K̄i et
l’intersection de deux triangles Ki et Kj avec i 6= j est soit vide, soit un point, soit une arête
entière. Les triangles (Ki )i=1,··· ,nk sont appelés les mailles ou les éléments du maillage. Pour
un triangle Ki , son diamètre est diam(Ki ) = max kx − yk. Il est appelé aussi taille de la
(x,y)∈Ki
maille Ki . La taille du maillage est définie comme h = max diam(Ki ).
Ω
On définit Pk l’ensemble des polynômes en 2D de degré inférieur ou égal à k ∈ N par :
X
Pk = {p(x, y) = aij xi y j , 0 ≤ i + j ≤ k, aij ∈ R} (2.1.14)
i,j≥0

On introduit alors l’espace de dimension finie Vh

Vh = {vh ∈ C(Ω̄); vh |K ∈ Pk , ∀K ∈ Th } (2.1.15)

La dimension de Vh est égale au nombre de noeuds de discrétisation du maillage. Ces noeuds

de discrétisation dépendent de k, degré des polynômes Pk choisi. Ils sont localisés sur les
sommets des triangles et éventuellement en d’autres points particuliers de ces triangles.
On définit alors Vh0 :
Vh0 = {vh ∈ Vh , vh = 0 sur ΓD } (2.1.16)
où la dimension de Vh0 égale à la différence entre la dimension de Vh et le nombre de noeuds
sur la frontière ΓD .
Ces deux espaces sont des sous-espaces de Sobolev de H 1 (Ω). Les fonctions de forme Ni
de Vh0 sont des fonctions "chapeau" qui prennent la valeur 1 en un noeud Si du maillage et
qui sont nulles en Sj avec i 6= j.

Ni (Sj ) = δij , 1 ≤ i, j ≤ n0h (2.1.17)

où n0h = nh − nD D
h avec nh le nombre de noeuds du maillage et nh les noeuds sur les bords ΓD .
2.1. Préliminaires 29

La base N = {N1 , · · · , Nn0 } associée aux noeuds du maillage est appelée base nodale.
h
On décompose uh sur cette base :

nh0
X
uh = Ni ui (2.1.18)
i=1

En prenant vh = Ni , ∀i ∈ {1, · · · , n0h }, le problème (2.1.13) s’écrit :

0
nh 0
nh
X X
a( N j u j , Ni ) = uj a(Nj , Ni ) = l(Ni ) ∀i = 1, · · · , nh . (2.1.19)
j=1 j=1

On peut alors écrire (2.1.19) sous la forme du système linéaire suivant :

0
Trouver U = (u1 , u2 , · · · , un0 )T ∈ Rnh tel que
(
h
(2.1.20)
KU = F

avec Kij = a(Nj , Ni ) et Fi = l(Ni ).

Après résolution du système linéaire (2.1.20), la solution (uh , σh ) est calculée par :

uh = N T U
(
(2.1.21)
σh = Cε(uh )

Remarque : dans X
le cas d’éléments de type quadrangle, l’espace d’approximation se base
sur Qk = {q(x, y) = aij xi y j , i ≤ k, j ≤ k, aij ∈ R}.
i,j≥0

Propriétés de la solution approchée par la méthode des éléments finis :

Pour la méthode des éléments finis de Lagrange que nous avons introduite ici, il s’agit de
chercher le déplacement dans un espace de dimension finie Vh0 . La solution (uh , σh ) obtenue
par la méthode des éléments finis est une approximation de la solution exacte (u, σ).
Le champ uh est dans H 1 . Il satisfait les conditions (2.1.7). Le champ de déplacement uh
est donc cinématiquement admissible.
1
Le champ σh = C( (graduh + gradT uh )) est dans L2 . Rien ne garantit que le vecteur σh n
2
soit continu le long des frontières des éléments. Le champ de contraintes σh n’est donc en
général pas statistiquement admissible. Il vérifie seulement l’équation équilibre (2.1.8) au
sens faible (intégration).
Le couple (uh , σh ) vérifie par contre la relation de comportement.

2.1.4 Erreur de discrétisation

La méthode des éléments finis est une des méthodes d’approximation les plus utilisées
dans le domaine de la mécanique des solides. Des outils ont été développés afin de permettre
d’estimer ses erreurs d’approximation et de les contrôler.
30 Chapitre 2. Raffinement de maillage en mécanique des solides

Notion d’erreur :
Quand on choisit d’utiliser une méthode numérique pour résoudre des problèmes physiques,
il existe différentes origines à l’erreur. Comme cela est montré dans [Ladevèze & Pelle 2005],
les principales sources d’erreur sont :
— Erreur de modèle : erreur de modélisation mathématique des problèmes physiques.
— Erreur de discrétisation : erreur entre la solution approchée et la solution exacte du
modèle.
— Erreur de résolution : erreur entre la solution approchée (définie par (2.1.20), par
exemple) et celle obtenue numériquement par une méthode de résolution.
Dans ce mémoire, on ne s’intéresse qu’à l’erreur de discrétisation et l’erreur de résolution.

Erreur de discrétisation :
On définit l’erreur de discrétisation en terme de déplacement comme

euh := u − uh (2.1.22)

On définit l’erreur de discrétisation en terme de contraintes :

eσh := σ − σh (2.1.23)

Normes pour mesurer l’erreur de discrétisation :

Nous introduisons différentes normes pouvant être utilisées pour mesurer les erreurs (2.1.22)
ou (2.1.23) :
Z 1/2
2
Norme L : kvkL2 := v · vdω
Ω
Z 1/2
1
Norme H : kvkH 1 := (v · v + ∇v : ∇v)dω
Ω
Z 1/2
Semi-norme H 1 : |v|H 1 := ∇v : ∇vdω
Ω

En mécanique des solides, nous utilisons souvent la "norme énergie" d’un champ de déplace-
ment qui est en fait une semi-norme équivalente à une norme lorsque le problème possède
des conditions aux limites de Dirichlet :
Z 1/2
Norme énergie : kvkE := σ(v) : ε(v)dω
Ω

Estimateur de l’erreur de discrétisation :

Il existe deux types d’estimateur d’erreur :
— Estimateur d’erreur a priori
L’estimateur d’erreur a priori fournit des informations utiles concernant le caractère
asymptotique de la solution approchée à partir de la géométrie du domaine, du
maillage utilisé, des fonctions d’interpolation de l’espace d’éléments finis, de la so-
lution exacte, etc. On pourra se reporter ici à l’ouvrage de Ciarlet [Ciarlet 1978].
2.2. Raffinement de maillage 31

Pour la méthode des éléments finis de Lagrange, l’estimateur a priori donne :

Si u ∈ H s (Ω), s > 1, il existe C > 0 tel que

∀h, ku − uh kH 1 (Ω) ≤ Chl |u|H l+1 (Ω) (2.1.24)

ou en norme L2 , ∀h, ku − uh kL2 (Ω) ≤ Chl+1 |u|H l+1 (Ω) (2.1.25)

avec l = min(s−1, k), k l’ordre d’interpolation de Lagrange utilisé. L’entier l est appelé
ordre de convergence de la méthode. Ainsi, on en déduit que :

lim ku − uh kH 1 (Ω) = 0 (2.1.26)

h→0

— Estimateur d’erreur a posteriori

L’estimateur d’erreur a posteriori ne dépend pas de la solution exacte. Il dépend des
données du problème, de la solution approchée obtenue, etc. Fournissant des infor-
mations sur la qualité de la solution approchée, il est très utilisé dans les techniques
de raffinement automatique de maillage. Ce type d’estimateur d’erreur englobe donc
l’erreur de discrétisation et l’erreur de résolution. Cependant, si le solveur utilisé est
suffisamment précis, l’erreur de résolution est en générale négligeable devant l’erreur
de discrétisation. Nous détaillons les différentes familles d’estimateurs d’erreur a pos-
teriori dans la section 2.3.

2.2 Raffinement de maillage

Dans le cadre de simulations industrielles, il est fréquent d’avoir des besoins en précision
différents selon les endroits du domaine de calcul. Dans notre cas d’étude, à savoir la modéli-
sation de l’interaction pastille-gaine, nous avons besoin d’une solution très précise autour des
zones de contact. La singularité de la solution dans ces zones oblige le maillage à être suffi-
samment fin. L’utilisation d’un maillage uniforme très fin sur tout domaine est très coûteuse
en terme de temps de calcul et d’espace mémoire. Le raffinement local de maillage fournit un
outil pour réaliser des simulations numériques de ce type de problème.
Dans cette section, nous présentons les principales méthodes de raffinement de maillage
s’appliquant à la méthode des éléments finis. Nous pouvons séparer les méthodes de raffine-
ment en deux familles : les méthodes adaptatives qui ont pour but de modifier le maillage
initial afin d’obtenir un maillage global optimal (taille d’élément, degré d’élément), et les
méthodes multi-grilles locales qui combinent (via un processus itératif) les solutions locales
obtenues sur plusieurs niveaux de maillages générés à partir du maillage initial.

2.2.1 Méthodes adaptatives

Les méthodes adaptatives reposent sur la modification locale du maillage initial. Nous
présentons les quatres types de méthodes adaptatives se différenciant par la manière dont le
maillage est modifié :
32 Chapitre 2. Raffinement de maillage en mécanique des solides

• Les méthodes h-adaptatives consistent à raffiner localement des éléments initiaux du

maillage grossier en éléments de taille de maille plus petite pour un degré de fonction
d’interpolation fixé (voir figure 2.1).
Un élément sera divisé en plusieurs éléments plus petits si l’erreur locale est trop im-
portante. Le rapport de raffinement peut varier selon les éléments. En utilisant un
raffinement par calcul de tailles d’éléments, le maillage obtenu est en général op-
timal immédiatement si l’indicateur d’erreur choisi est fiable. Différentes méthodes
de calcul de tailles optimales d’éléments ont été proposées [Zienkiewicz & Zhu 1987,
Ladevèze et al. 1991, Oñate & Bugeda 1993, Li et al. 1995]. Une comparaison des per-
formances de ces méthodes, étudiée par [Ehlers et al. 2002], montre que la méthode
de [Ladevèze et al. 1991] conduit au maillage le plus optimal.

F I G U R E 2.1 – Exemple de h-raffinement d’un maillage de quadrangles

Une autre stratégie consiste à raffiner récursivement les mailles avec un rapport de
raffinement constant (voir la figure 2.2).

F I G U R E 2.2 – Exemple de raffinement h-adptatif récursif

Les méthodes h-adaptatives permettent d’améliorer la solution quand le problème

présente des singularités [Babuska & Szabo 1982]. Ces méthodes sont les mé-
thodes de raffinement les plus utilisées (par exemple, [Demkowicz et al. 1985,
Strouboulis & Haque 1992, Belytschko & Tabbara 1993, Fish & Markolefas 1994,
Babuška et al. 1994, Rachowicz 1997, Han & Wriggers 2000, Cai & Li 2010,
Bao et al. 2012]).
Quelle que soit la stratégie choisie, les méthodes h-adaptatives conduisent générale-
ment à des maillages non-conformes, surtout pour des quadrangles. Il faut alors soit
définir des relations supplémentaires sur les noeuds non-conformes afin de permettre
de résoudre le problème défini sur ce maillage par la méthodes des éléments finis, soit
trouver des procédures permettant de rendre le maillage conforme. La première solu-
tion nécessite de disposer d’un solveur éléments finis adapté. Plusieurs procédures ont
ainsi été développées afin de disposer d’un maillage raffiné conforme :
2.2. Raffinement de maillage 33

— Dans le cas d’éléments quadrangulaires, Strouboulis et Haque

[Strouboulis & Haque 1992] ont proposé de définir une zone de transition
sous forme d’un mélange entre des éléments triangulaires et quadrangulaires
(voir figure 2.3). Mais il faut alors que le solveur éléments finis puisse traiter des
maillages avec différentes types d’éléments.

h-méthode
local conforme

Zone d'intérêt

: noeud d'interpolation

F I G U R E 2.3 – h-raffinement conforme avec différent type d’éléments

— Pour les éléments triangles, plusieurs algorithmes ont été proposés afin de rendre
les maillages conformes. Par exemple, l’algorithme basé sur la bissection du bord
le plus long proposé par Rivara [Rivara 1984], l’algorithme récursif Newest Vertex
Bisection proposé par Mitchell [Mitchell 1991] (voir exemple de figure 2.4a), l’al-
gorithme de raffinement d’un triangle en triangles simulaires proposé par Bank
[Bank et al. 1983] (voir exemple de figure 2.4b), ...

(a) (b)

F I G U R E 2.4 – h-raffinement conforme d’un maillage triangulaire

• Les méthodes p-adaptatives cherchent à augmenter localement le degré des fonctions

d’interpolation tandis que la géométrie du maillage reste inchangée (voir figure 2.5).

p-méthode

Zone d'intérêt

: noeud d'interpolation

F I G U R E 2.5 – Méthode p-raffinement

La vitesse de convergence des méthodes p-adaptatives est souvent plus élevée que
celle des méthodes h-raffinement et peut être exponentielle [Babuska et al. 1981,
34 Chapitre 2. Raffinement de maillage en mécanique des solides

Babuska & Szabo 1982]. Ce type de méthodes est également intéressant en terme de
temps de génération de maillage. Elles sont aussi beaucoup étudiées dans la litté-
rature (par exemple, [Babuška & Suri 1987a, Düster & Rank 2001, Barros et al. 2004,
Kubatko et al. 2009]). La plupart des implémentations de ce type de méthodes utilisent
des fonctions de forme hiérarchiques cf. la figure 2.6 (voir [Zienkiewicz et al. 1983]
pour plus détails) permettant simplement d’obtenir le degré d’approximation désiré
en ajoutant de nouvelles fonctions de formes aux fonctions de formes existantes. Les
fonctions de forme obtenues sont de plus généralement conformes.

1 1

0 0
noeud noeud
Fonctions de forme d'ordre 1 Fonctions de forme d'ordre 2

F I G U R E 2.6 – Exemple des fonctions de forme hiérarchiques en 1D

Cependant, les méthodes p-adaptatives ne permettent pas d’améliorer l’approxima-

tion de la géométrie du problème car elles s’appuient toujours sur le maillage ini-
tial du domaine. Elles ne permettent pas forcément de mieux approximer des sin-
gularités de solution. Elles sont par ailleurs peu utilisées dans le domaine industriel
car il faut disposer d’un solveur sachant gérer des éléments d’ordres différents et les
codes industriels utilisent rarement des éléments dont l’ordre est supérieur à 2 (par
exemple,[Ladevèze & Pelle 2005, CAST3M 2016]).
• Les méthodes r-adaptatives consistent à déplacer des noeuds du maillage vers les zones
d’intérêt (voir figure 2.7).

r-méthode

Zone d'intérêt

: noeud d'interpolation

F I G U R E 2.7 – Méthodes r-raffinement

Les méthodes r-adaptatives sont intéressantes pour des problèmes transitoires où les élé-
ments peuvent suivre des phénomènes en évolution. Ce type de méthodes est aussi beau-
coup étudié (par exemple [Ghosh & Manna 1993, Cao et al. 2001, Beckett et al. 2001,
Huang & Russell 2010]). Différents algorithmes en déplacement ont été développés,
le lecteur pourra notamment se référer à [Huang & Russell 2010]. Les méthodes r-
adaptatives gardent la conformité du maillage et sont simples à mettre en oeuvre. Ce-
2.2. Raffinement de maillage 35

pendant, l’amélioration de la solution est conditionnée par le nombre initial de noeuds

du maillage.
• Les méthodes s-adaptatives, proposées par [Fish 1992], superposent des maillages plus
fins sur le maillage grossier dans les zones d’intérêt afin de former un maillage compo-
site (voir la figure 2.8). Sur les bords "fictifs" des maillages fins ∂ΩL \∂Ω, des conditions
aux limites de Dirichlet homogènes sont imposées. Si les maillages fins ont des bords
en communs avec les bords de Dirichlet du maillage grossier ∂ΩL ∩ ΓD 6= ∅, des condi-
tions aux limites de Dirichlet homogènes y sont également imposées. Les degrés de
liberté des maillages fins qui sont communs au maillage grossier et aux maillages fins
sont supprimés. Après la résolution d’un système composite, la solution est la somme
des solutions de tous les niveaux. Ces méthodes ont été appliquées à des problèmes
linéaires et nonlinéaires (par exemple,[Fish & Markolefas 1993, Fish & Guttal 1996,
Yue & Robbins 2005, Yue & Robbins 2007]). Cependant, elles augmentent rapidement
la dimension de la matrice du systèmes linéaire à inverser.
Maillage grossier

Maillage fin
Ω
ΩL
ΩL

Maillage composite

Ω
ΩL

F I G U R E 2.8 – Exemple simple de superposition des maillages [Yue & Robbins 2005]

La combinaison de différentes méthodes de raffinement permet de cumuler

leurs avantages. Des méthodes hp-adaptatives (par exemple,[Guo & Babuška 1986,
Babuška & Suri 1987b, Demkowicz et al. 1989, Rachowicz et al. 1989, Oden & Demkowicz 1991,
Rachowicz & Zdunek 2005, Düster et al. 2007, Solin et al. 2010, Lebrun-Grandié et al. 2011])
et des méthodes hr-adaptatives (par exemple, [Sun & Zamani 1992, Edwards et al. 1993,
Askes & Rodriguez-Ferran 2001]) ont notamment été développées.

2.2.2 Méthodes multi-grilles locales

Les méthodes multi-grilles standards (voir [Fedorenko 1961, Brandt 1977, Hackbusch 1985])
permettent d’accélérer la convergence de la résolution d’un système contenant beaucoup de
36 Chapitre 2. Raffinement de maillage en mécanique des solides

degrés de liberté. Elles se basent sur le fait que les méthodes itératives éliminent rapidement
les hautes fréquences de l’erreur et que les parts de basses fréquences de l’erreur peuvent être
bien calculées sur des maillages grossiers. A partir d’un maillage initial fin, plusieurs niveaux
de grilles grossières sont définis sur tout le domaine (voir la figure 2.9). Des opérateurs de
prolongement et de restriction relient les différentes grilles jusqu’à la convergence de la solu-
tion sur le maillage initial fin. L’algorithme multi-grilles est une généralisation de l’algorithme
bi-grilles. L’algorithme bi-grilles consiste à effectuer itérativement jusqu’à convergence de
la solution fine quelques itérations de lissage sur le système à résoudre sur le maillage fin,
puis à restreindre sur le maillage grossier les résidus fins obtenus et utiliser cette restriction
comme correction du second membre du problème grossier. La solution exacte du problème
grossier est quant à elle ensuite prolongée sur le maillage fin pour corriger la solution fine.
L’opérateur de prolongement doit être la transposée de l’opérateur de restriction pour assurer
la convergence des méthodes multi-grilles.
Il est à noter que ce type de méthodes a été très étudié dans la litté-
rature (par exemple,[Lubrecht et al. 1987, Lebon 1995, Boersma & Wriggers 1997,
Wriggers & Boersma 1998, Lebon et al. 2007]) et a été notament utilisé pour la résolu-
tion de problèmes de contact avec ou sans frottement, qui feront l’objet du Chapitre 4.

Niveau 2
Niveau 0

Niveau 1
Niveau 1

Niveau 2 Niveau 0

F I G U R E 2.9 – Exemple de 3 niveaux F I G U R E 2.10 – Exemple de 3 niveaux de

de grilles du processus multi-grilles grilles du processus multi-grilles local

Toutes les méthodes multi-grilles locales sont issus des travaux réputés de Brandt
[Brandt 1977] dans le cadre du développement des méthodes MLAT (Multi-Level Adaptive
Techniques). Les méthodes multi-grilles locales elles aussi consistent à résoudre un problème
sur plusieurs niveaux de grilles (voir la figure 2.10). Cependant, seul le maillage le plus gros-
sier recouvre le domaine initial de calcul. A partir de ce maillage grossier défini sur l’ensemble
du domaine, ces méthodes consistent à ajouter récursivement des niveaux locaux avec pas de
mailles de plus en plus fins dans des sous-domaines où la solution n’est pas suffisamment pré-
cise. Un processus multi-grille (voir la figure 2.11 comme exemple) relie les différents niveaux
2.2. Raffinement de maillage 37

de grilles jusqu’à convergence de la solution grossière. L’opérateur de prolongement consiste

à définir (souvent par interpolation) des conditions aux limites de Dirichlet sur les bords fictifs
des maillages fins à partir des solutions grossières. L’opérateur de restriction permet quant
à lui de corriger les solutions sur les maillages grossiers à partir de résidus obtenus via les
solutions fines. Les opérateurs de restriction varient selon la méthode multi-grilles locales
choisie.
Les méthodes multi-grilles locales sont des méthodes très génériques. Le maillage, le
solveur, le modèle et le rapport de raffinement peuvent changer entre les niveaux.
Grille fine

Grille grossière

: Résolution

F I G U R E 2.11 – Exemple de processus multi-grilles local pour 4 niveaux de grille : ∧-cycles

Dans la littérature, il existe plusieurs méthodes de raffinement multi-grilles local. Nous

présentons ici celles souvent utilisées avec la méthode des éléments finis.
• La méthode LDC (Local Defect Correction) a été proposée par Hackbusch
[Hackbusch 1984]. L’opérateur de restriction de cette méthode consiste à restreindre
la solution fine sur le maillage grossier afin de calculer un résidu grossier associé
sur la zone de recouvrement. Le problème grossier est alors résolu en ajoutant
au second membre le résidu obtenu sur les zones de recouvrement. Ferket et al.
[Ferket 1996, Ferket & Reusken 1996b, Ferket & Reusken 1996a] ont complété l’ana-
lyse mathématique de Hackbusch dans le cadre d’une discrétisation par différences
finies. La combinaison de la méthode LDC avec la méthode des volumes finis a
été étudiée par Anthonissen [Anthonissen et al. 1998]. L’analyse de la méthode
LDC combinée avec la méthode des éléments finis a été effectuée par Wappler
[Wappler 1999]. Cette méthode, beaucoup étudiée en mécanique des fluides (par
exemple, [Belliard & Grandotto 2003, Anthonissen et al. 2005, Minero et al. 2006,
Kramer et al. 2009]) a été récemment appliquée avec succès à la mécanique des solides
[Barbié et al. 2014, Barbié et al. 2015].
• Les méthodes "Full MultiGrid" (FMG) locales [Cavin et al. 2005, Watremetz et al. 2007,
Biotteau 2010] cherchent à générer et relier les sous-niveaux locaux par un processus
de type FMG [Brandt 1994] (voir figure 2.12). Le processus FMG consiste à accélérer la
vitesse de convergence de la méthode multi-grilles en faisant démarrer l’algorithme du
maillage le plus grossier et en effectuant une résolution bi-grilles. Les niveaux fins sont
ajoutés successivement après convergence des cycles multi-grilles précédents. Ainsi
dans le processus FMG local, les sous-niveaux les plus fins ne sont ajoutés que si la
38 Chapitre 2. Raffinement de maillage en mécanique des solides

solution corrigée n’est pas suffisamment précise. L’opérateur de restriction est identique
à celui de la méthode LDC. L’opérateur de prolongement consiste en plus à définir
une estimation initiale de la solution fine à partir de l’interpolation de la solution
grossière. La méthode FMG locale ressemble donc beaucoup à la méthode LDC. La
seule différence réside dans le fait que les sous-niveaux FMG locaux ne sont ajoutés
qu’après convergence des solutions sur les niveaux précédents.

Grille fine

Grille grossière

: Résolution

F I G U R E 2.12 – Représentation du cycle FMG pour 4 niveaux de grilles

• La méthode FAC (Fast Adaptive Composite), proposée par McCormick [McCormick 1984,
McCormick & Thomas 1986], introduit un problème composite intermédiaire entre
deux niveaux de grille. Le résidu calculé sur le maillage composite est utilisé pour
corriger à la fois le problème grossier et le problème fin. Cette méthode est souvent
utilisée avec une discrétisation par différences finies. L’analyse de convergence de cette
méthode a été étudiée dans [McCormick 1984, McCormick & Thomas 1986]. Les résul-
tats numériques (par exemple, [Hart et al. 1986, Thomas et al. 1986, Cai et al. 1995,
Ritter et al. 2010]) montrent l’efficacité de cette méthode qui est par ailleurs facilement
parallèlisable.
Il est à noter que la méthode FAC a été aussi appliquée à des problèmes de contact (voir
[Grego 1995]).

2.3 Estimateurs d’erreur a posteriori

2.3.1 Généralités
Les estimateurs d’erreur a posteriori sont des outils indispensables aux techniques de raffi-
nement automatique de maillage. Ils fournissent des informations sur la qualité de la solution
approchée qui permettent de détecter les zones d’intérêt dans le processus de raffinement
automatique.
Dans cette section, nous présentons les trois principales familles d’estimateurs d’erreur a
posteriori :
• Les estimateurs d’erreur basés sur la vérification de la loi de comportement, introduits
par Ladevèze [Ladevèze 1975].
2.3. Estimateurs d’erreur a posteriori 39

• Les estimateurs d’erreur basés sur le résidu de l’équation d’équilibre, introduits par
Babuška et Rheinboldt [Babuška & Rheinboldt 1978b].
• Les estimateurs d’erreur basés sur le lissage des contraintes, introduits par Zienkiewicz
et Zhu [Zienkiewicz & Zhu 1987].

Tous les estimateurs précédents reposent sur le fait que comme indiqué dans le paragraphe
2.1.3 en page 29, la solution approchée par la méthode des éléments finis n’est pas statique-
ment admissible.
Ladevèze [Ladevèze 1975] a proposé un estimateur basé sur l’erreur en relation de com-
portement. Il cherche alors à construire une solution statiquement admissible à partir de
la solution cinématiquement admissible obtenue par la méthode des éléments finis. Ce-
pendant, il n’y a pas de raison que cette nouvelle solution vérifie la loi de comporte-
ment. L’erreur a posteriori est alors définie comme l’erreur sur la relation de comporte-
ment de cette solution statiquement et cinématiquement admissible. Des résultats d’utili-
sation de cette méthode en 2D élastique pour les éléments finis P1 ont été publiés dans
[Ladevèze & Leguillon 1983, Ladevèze et al. 1986]. Dans [Ladevèze et al. 1991], Ladevèze a
étendu cette méthode à d’autres types d’éléments finis pour des problèmes élastiques en 2D.
Le point de départ du développement de la seconde famille d’estimateurs d’erreur, les esti-
mateurs d’erreur basés sur le résidu de l’équation d’équilibre, est que la solution éléments
finis ne vérifie pas localement l’équation d’équilibre (2.1.1). En se basant sur le calcul du
résidu de l’équation d’équilibre du problème mécanique linéaire 1D, Babuška et Rheinholdt
[Babuška & Rheinboldt 1978b] ont développé des techniques d’estimation d’erreur a poste-
riori qui permet d’approximer l’erreur sur chaque élément fini K en norme énergie.
Plus récemment, Zienkiewicz et Zhu [Zienkiewicz & Zhu 1987] ont développé une tech-
nique simple afin d’estimer l’erreur de discrétisation. Cette technique est basée sur le fait
que les contraintes obtenues par la méthodes des éléments finis sont discontinues entre
les éléments lorsque on utilise des éléments finis de Lagrange car uh est dans l’espace
H 1 tandis que σh est dans L2 . Les auteurs ont proposé de lisser les contraintes des élé-
ments finis pour obtenir des contraintes plus régulières et d’estimer l’erreur à partir de
la différence entre ces contraintes lissées et les contraintes obtenues par la méthode des
éléments finis. Ainsworth et al. [Ainsworth et al. 1989] ont montré que la solution lissée
converge vers la solution exacte si la solution exacte est suffisamment régulière. Plus tard,
la méthode SPR (Superconvergent Patch Recovery) a été proposée par Zienkiewicz et Zhu
[Zienkiewicz & Zhu 1992a, Zienkiewicz & Zhu 1992b]. Cette méthode consiste à construire
sur des patchs d’éléments le champ de contraintes lissées à partir de points dits "superconver-
gents".
La définition mathématique de ces estimateurs va être donnée dans les sections
2.3.3 à 2.3.5. Cependant, pour plus de détails, le lecteur peut aussi se reporter à
[Ladevèze & Pelle 2005, Ainsworth & Oden 1997, Gratsch & Bathe 2005, Verfürth 1999,
Ramm et al. 2003]. La relation entre les différents estimateurs a été étudiée par Zhu
[Zhu 1997]. Parmi les différentes familles d’estimateur, les estimateurs basés sur le lissage
des contraintes sont les plus utilisés.
40 Chapitre 2. Raffinement de maillage en mécanique des solides

2.3.2 Critères de comparaison des estimateurs

Afin de comparer les performances des différentes familles d’estimateurs a posteriori, des
critères de comparaison doivent être définis. On en introduit plusieurs dans ce qui suit.
— Indice d’efficacité :
Soit ηK la norme d’un estimateur d’erreur local sur un élément K, on peut calculer la
norme globale η de l’estimateur d’erreur avec
 1/2
X 
2
η= ηK
 
K∈Th

L’indice d’efficacité est le rapport entre la norme globale η de l’estimateur d’erreur et

la norme de l’erreur de discrétisation keh k :
η
γ=
keh k
où l’erreur de discrétisation globale eh est définie comme la différence entre la solution
analytique et la solution numérique (voir (2.1.22) ou (2.1.23)), k·k désigne une norme
(choix arbitraire). En général, on ne connaît pas la solution analytique. On utilise donc
une solution approchée sur un maillage très fin pour calculer l’erreur de discrétisation.
Un bon estimateur doit avoir un indice d’efficacité proche de 1, même pour un maillage
assez grossier, afin d’être utilisable.
— Fiabilité :
On dit qu’un estimateur est fiable si l’erreur estimée associée tend vers zéro quand la
taille du maillage tend vers zéro. Idéalement, l’erreur estimée et l’erreur de discréti-
sation doivent tendre vers zéro simultanément. De plus, l’erreur estimée doit garantir
des bornes supérieure et inférieure de l’erreur de discrétisation.
— Simplicité d’implantation et coût d’utilisation
Un bon estimateur doit être simple à implémenter et rapide à utiliser afin de ne pas
polluer le temps de calcul global de la méthode de résolution.

2.3.3 Estimateur basé sur le résidu d’équilibre

Il existe deux types d’estimateurs dans la famille des estimateurs basés sur le résidu
d’équilibre : des estimateurs dits explicites et des estimateurs dits implicites. Les estimateurs
d’erreur basés directement sur la solution obtenue par la méthode des éléments finis sont
appelés estimateurs explicites. Au contraire, les estimateurs d’erreur implicites ont besoin
de résoudre des problèmes locaux auxiliaires. En conséquence, les estimateurs explicites
nécessitent un coût de calcul moindre mais sont souvent moins précis que les estimateurs
implicites.

[Link] Estimateurs d’erreur explicites

L’estimateur d’erreur explicite est calculé directement à partir des résidus élémen-
taires de l’équation d’équilibre et des sauts des contraintes aux interfaces des éléments.
2.3. Estimateurs d’erreur a posteriori 41

Ces techniques d’estimation ont été proposées initialement par Babuška et Rheinboldt
[Babuška & Rheinboldt 1978a, Babuška & Rheinboldt 1981].
Par définition (cf. éq. (2.1.22)), euh satisfait :

a(euh , v) = a(u, v) − a(uh , v) ∀v ∈ V 0 (2.3.1)

D’après (2.1.11), on a :

a(euh , v) = l(v) − a(uh , v) ∀v ∈ V 0 (2.3.2)

On décompose l’équation (2.3.2) sur chaque élément :

a(euh , v) = l(v) − a(uh , v)

X Z Z Z 
= f · v dΩ + FN · v dΓ − Cε(uh ) : ε(v) dΩ ∀v ∈ V 0
K∈Th K ∂K∩ΓN K

(2.3.3)

avec ∂K la frontière de l’élément K.

On note :
∂Th = ∪K∈Th ∂K.
En intégrant le dernier terme de l’équation (2.3.3) par parties, on obtient :
X Z X Z
u
a(eh , v) = R · v dΩ + J · vdΓ ∀v ∈ V 0 (2.3.4)
K∈Th K ∂K∈∂Th ∂K

où R est le résidu élémentaire :

R = f + divσh dans K (2.3.5)

et J le saut des contraintes à travers ∂K :


′ ′
−(σh n + σh n ) si ∂K ∈
/Γ


J= FN − σ h n si ∂K ∈ ΓN (2.3.6)


0 si ∂K ∈ ΓD


où à l’intérieur du maillage, la frontière ∂K sépare les éléments K et K ′ (voir figure 2.13).

La condition d’orthogonalité de Galerkin donne :

a(euh , vh ) = 0 ∀vh ∈ Vh0 ⊂ V 0 (2.3.7)

On a ainsi

a(euh , v) = a(euh , v) − a(euh , vh ) = a(euh , v − vh ) ∀vh ∈ Vh0 , ∀v ∈ V 0 (2.3.8)

En introduisant l’interpolation de Clément [Clément 1975] :

Soit χh un opérateur d’interpolation linéaire de V 0 → Vh0 : v 7→ χh v (2.3.9)

42 Chapitre 2. Raffinement de maillage en mécanique des solides

F I G U R E 2.13 – Interface entre deux éléments

on a :
X Z X Z
a(euh , v) = a(euh , v − χh v) = R (v − χh v) dΩ + J(v − χh v)dΓ ∀v ∈ V 0
K∈Th K ∂K∈∂Th ∂K
(2.3.10)
On applique ensuite l’inégalité de Cauchy-Schwarz :
X X
a(euh , v) ≤ kRkL2 (K) kv − χh vkL2 (K) + kJkL2 (∂K) kv − χh vkL2 (∂K) (2.3.11)
K∈Th ∂K∈∂Th

Selon les résultats du théorème d’interpolation de Clément [Clément 1975], nous avons :

kv − χh vkL2 (K) ≤ c01 hK kvkH 1 (K̃) (2.3.12)

p
kv − χh vkL2 (∂K) ≤ c02 hK kvkH 1 (K̃) (2.3.13)
avec hK le diamètre de l’élément K, c01 , c02 des constantes qui dépendent de la forme des
éléments, K̃ est l’ensemble des éléments qui ont une arête commune avec K.
Par l’inégalité de Hölder, on obtient :
X X p
a(euh , v) ≤ c11 kRkL2 (K) hK kvkH 1 (K̃) + c22 kJkL2 (∂K) hK kvkH 1 (K̃)
K∈Th ∂K∈∂Th
 1/2
X X 
≤ c33 max{c11 , c22 }kvkH 1 (Ω) h2K kRk2L2 (K) + hK kJk2L2 (∂K)
 
K∈Th ∂K∈∂Th
(2.3.14)
avec c11 = max c01 et c22 = max c02 .
K∈Th K∈Th
En remplaçant v par eh et grâce à la définition de la norme énergie, on a :
 1/2
X X 
keuh k2E = a(euh , euh ) ≤ c33 max{c11 , c22 }keuh kH 1 (Ω) h2K kRk2L2 (K) + hK kJk2L2 (∂K)
 
K∈Th ∂K∈∂Th
(2.3.15)
L’équivalence de norme (c0 kvkH 1 (Ω) ≤ kvkE(Ω) ≤ kvkH 1 (Ω) avec constante c0 > 0) donne :
 
X X 
keuh k2E ≤ c h2K kRk2L2 (K) + hK kJk2L2 (∂K) (2.3.16)
 
K∈Th ∂K∈∂Th
2.3. Estimateurs d’erreur a posteriori 43

1
avec c = ( c33 max{c11 , c22 })2 .
c0
On a obtenu un estimateur d’erreur globale basé sur le résidu de l’équation d’équilibre :
 
X X 
η2 = c h2K kRk2L2 (K) + hK kJk2L2 (∂K) (2.3.17)
 
K∈Th ∂K∈∂Th

L’indicateur d’erreur locale pour chaque élément ηK peut être calculé par :
n o
2
ηK = ch2K kRk2L2 (K) + chK kJk2L2 (∂K) (2.3.18)

et X
η2 = 2
ηK (2.3.19)
K∈Th

Ce type d’estimateur est mathématiquement rigoureux. L’estimateur garantit la borne supé-

rieure de l’erreur. Ces estimateurs sont relativement simples à mettre en oeuvre.
Par contre, l’indice d’efficacité est relié à la constante c qui est difficile à calculer. Les
estimateurs sont donc généralement pessimistes car c est souvent dictée par des critères de
sécurité [Gratsch & Bathe 2005]. Ceci explique que cet estimateur sous la forme (2.3.18) est
peu utilisé dans la littérature.
Le cas 2D a été étudié dans les travaux de Babuška et Miller [Babuška & Miller 1987] et
Kelly et al. [Kelly et al. 1983]. Kelly et al. [Kelly et al. 1983] ont donné une forme explicite
de la constante c pour les problèmes en 1D et en 2D. Pour les éléments Pk en 2D, l’estimateur
s’écrit :  2 
2 hK 2 hK 2
ηK = kRkL2 (K) + kJkL2 (∂K) (2.3.20)
24kρ 24kρ
où k est le degré d’interpolation de l’élément fini et ρ est une constante dépendant du ma-
tériau. Pour le cas des contraintes planes, Kelly et al. [Kelly et al. 1983] ont montré que
E
m= avec E le module d’Young, ν le coefficient de Poisson. Cette formule est souvent
1−ν
utilisée [Zienkiewicz et al. 2005]. Elle est implantée dans Code_Aster [Delmas 2011]. Cepen-
dant, on ne dispose pas d’évaluation de la constante pour le cas 3D qui nous intéresse pour la
modélisation 3D de l’IPG.
Verfürth [Verfürth 1996] a développé une technique pour obtenir les bornes de l’er-
reur. L’estimateur en norme infinie a été analysé par Nochetto [Nochetto 1995]. Dari et al.
[Dari et al. 1999] ont étendu ce type d’estimateur en norme infinie au cas 3D.

[Link] Estimateurs d’erreur implicites

Du fait des difficultés de calcul de la constante c pour les estimateurs explicites, des esti-
mateurs d’erreur implicites ont été développés. Par la résolution d’un problème auxiliaire, les
estimateurs implicites permettent d’obtenir des informations plus précises que les estimateurs
explicites. Il existe deux types de méthodes : la méthode du résidu élémentaire et la méthode
du résidu par patch.
44 Chapitre 2. Raffinement de maillage en mécanique des solides

Méthode du résidu élémentaire :

La méthode du résidu élémentaire proposé par Bank and Weiser [Bank & Weiser 1985]
consiste à chercher l’erreur locale eK = euh|K = (u − uh )|K sur l’élément K qui satisfait
le problème variationnel suivant :
Z Z
a(eK , v)|K = R · v dΩ + (FN |K − σh|K n) · vdΓ
K ∂K∩ΓN
Z (2.3.21)
+ (σ|K n − σh|K n) · vdΓ ∀v ∈ VK
∂K\Γ

avec l’espace de fonctions test VK qui est défini par :

VK = {v ∈ H 1 (K) : v = 0 sur ∂K ∩ ΓD } (2.3.22)

Les deux premiers termes sont entièrement calculables à partir de la solution éléments
finis. Par contre, le troisième terme est inconnu car on ne connaît pas la solution exacte σ|K .
Celle-ci est approchée par :
1
σ|K n ≈ (σh|K + σh|K ′ )n (2.3.23)
2
′
où σh|K ′ est la contrainte dans l’élément K adjacent à l’élément K.
Le problème local approché devient :


 Trouver ζK ∈ VK tel que

 Z Z

a(ζK , v)|K = (fh|K + divσh|K ) · v dΩ + (Fn|K − σh|K n) · vdΓ

K ∂K∩ΓN (2.3.24)

1

 Z
+ (σh|K ′ n − σh|K n) · vdΓ ∀v ∈ VK


∂K\Γ 2


Après la résolution de (2.3.24) sur chaque élément K, l’erreur peut être estimée par :
X
keuh k2E ≈ 2
ηK avec ηK2
= kζK k2E (2.3.25)
K∈Ωh

Malheureusement, l’existence et l’unicité du problème (2.3.24) ne sont pas garanties à

cause des conditions aux limites de Neumann du problème local (quand ∂K ∩ ΓD = ∅). Cer-
taines méthodes sont proposées dans [Ainsworth & Oden 2011, Babuška & Strouboulis 2001]
pour pallier ce problème. Nous ne les détaillons pas ici.

Méthode du résidu par patch :

L’idée de la méthode de Babuška [Babuška & Rheinboldt 1978b, Babuška & Rheinboldt 1978a]
est de décomposer l’équation (2.3.2) en plusieurs problèmes locaux sur un patch d’éléments
avec des conditions aux limites de Dirichlet homogènes. Les autres étapes sont identiques
à celles de la méthode du résidu élémentaire. Plus récemment, d’autres méthodes ont été
proposés [Carstensen & Funken 1999, Prudhomme et al. 2004] qui sont plus flexibles par
rapport aux conditions aux limites. Prudomme et al. [Prudhomme et al. 2004] proposent
d’utiliser la partition de l’unité de la fonction de forme. Nous ne détaillons pas ici ce type de
méthode.
2.3. Estimateurs d’erreur a posteriori 45

En utilisant les méthodes des estimateurs implicites, on s’affranchit de l’étape de détermi-

nation de la valeur de la constante c. Par contre, il est obligatoire de résoudre des problèmes
locaux. Ces méthodes sont beaucoup plus coûteuses que des méthodes explicites mais elles
donnent des résultats plus précis. Ce type d’estimateur a été appliqué pour résoudre les
problème de Stokes par Bank et al. [Bank & Welfert 1991], Oden et al. [Oden et al. 1994].
Rachowicz et Zdunek [Rachowicz & Zdunek 2005] ont appliqué ce type d’estimateur avec la
méthode hp-adaptative sur les équations de Maxwell en 3D.

2.3.4 Estimateur basé sur la loi de comportement

Sachant que la solution obtenue par éléments finis est cinématiquement admissible mais
qu’elle n’est pas statiquement admissible (régularité du champ des contraintes), les estima-
teurs basés sur la loi de comportement consistent à construire une solution (uCA , σSA ) qui
est à la fois cinématiquement admissible (CA) et statiquement admissible (SA) et mais qui ne
respecte pas forcément la loi de comportement. Comme la solution en déplacement obtenue
avec la méthode des éléments finis est CA, on choisit uCA = uh .
La difficulté de ce type d’estimateur est donc de construire le champ de
contraintes statiquement admissible. De nombreuses techniques ont été déve-
loppées. Les trois principales techniques sont la technique d’équilibre par élé-
ment (EET) [Ladevèze & Leguillon 1983, Ladevèze et al. 1991, Ladevèze et al. 1986,
Ladevèze & Pelle 2005], la technique d’équilibrage par patch d’éléments SPET (ou
méthode flux-free) [Pares et al. 2006, Machiels et al. 2000, Carstensen & Funken 1999,
Gallimard 2009] et la technique d’équilibrage par éléments et patch d’éléments (ou mé-
thode hybride) [Ladevèze et al. 2010]. La comparaison des différentes techniques a été
publiée par Pled et al. [Pled et al. 2011]. Les auteurs ont montré que la méthode SPET fournit
l’estimateur le plus précis même si les trois techniques garantissent les bornes supérieures de
l’erreur en norme energie.
La méthode standard [Ladevèze & Leguillon 1983], aussi appelée la technique d’équili-
brage par élément EET, est basée sur la condition de prolongement :
Z
(σSA − σh ) gradNi dΩ = 0 ∀K ∈ Th , ∀i ∈ I (2.3.26)
K

où I = {1, ..., n0h } est l’ensemble des noeuds du maillage et Ni la fonction de forme associée
au noeud i . Cette condition de prolongement permet de prolonger σh en σSA .
Il y a deux étapes successives pour construire le champ statiquement admissible :
— construction du second membre élémentaire F̄K à partir de la solution éléments finis
pour chaque élément K.
— calcul des champs de contraintes admissibles par résolution des problèmes d’équilibre
locaux avec F̄K donnée.
Ce processus est maintenant décrit :
46 Chapitre 2. Raffinement de maillage en mécanique des solides

Construction du second membre élémentaire :

Les contraintes statiquement admissibles sont supposées respecter la condition :

σSA|K nK = αF̄K sur ∂K, ∀K ∈ Th (2.3.27)

avec α un entier relatif défini sur chaque face de l’élément et valant +1 ou −1 et F̄K le second
membre recherché.
Les fonctions αF̄K et σSA|K satisfont l’équation d’équilibre au sens faible :
Z Z Z
σSA gradNi dΩ = αF̄K · Ni dΓ + f · Ni dΩ ∀K ∈ Th , ∀i ∈ I (2.3.28)
K ∂K K

En utilisation la condition de prolongement (2.3.26), on a :

Z Z Z
σh gradNi dΩ = αF̄K · Ni dΓ + f · Ni dΩ ∀K ∈ Th , ∀i ∈ I (2.3.29)
K ∂K K

On obtient :
Z Z Z
αF̄K · Ni dΓ = σh gradNi dΩ − f · Ni dΩ ∀K ∈ Th , ∀i ∈ I (2.3.30)
∂K K K
Z
On a donc la valeur de αF̄K · Ni dΓ pour les frontières ∂K connectées au noeud i.
∂K

Résolution des problèmes élémentaires :

On cherche à résoudre les problèmes élémentaires :
(
− divσSA = f|K dans K
(2.3.31)
σSA n = αF̄K sur ∂K

On peut résoudre ce problème en utilisant la méthode des éléments finis en contraintes. Les
conditions aux limites sont alors des conditions de Dirichlet. On discrétise chaque élément de
la triangulation en sous éléments finis.
La solution (uCA , σSA ) est à la fois cinématiquement admissible et statiquement admissible,
mais elle n’a aucune raison de vérifier la loi de comportement.
On peut alors estimer l’erreur en relation de comportement par :

eRdC = σSA − Cε(uCA ) (2.3.32)

Les estimateurs basés sur le concept d’erreur en relation de comportement sont précis
avec des bornes garanties de l’erreur globale [Ladevèze & Pelle 2005]. Ces estimateurs ont
par exemple été appliqués à des modèles industriels : un modèle 2D de structure fissurée
[Pares et al. 2006, Ladevèze et al. 2010], un modèle 3D de plaque trouée en traction et un
modèle 3D de tête de réacteur [Pled et al. 2011] ... Par contre la construction de la solution
statiquement admissible est coûteuse.
Une méthode similaire mais basée sur le calcul direct des champs cinématiquement et
statiquement admissibles a aussi été étudiée par plusieurs auteurs [Coorevits et al. 2000,
Bellenger & Coorevits 2005, Kuss & Lebon 2011] pour les problèmes élastiques.
2.3. Estimateurs d’erreur a posteriori 47

2.3.5 Estimateur basé sur le lissage des contraintes

Les contraintes calculées par la méthode des éléments finis de Lagrange présentent des
discontinuités sur les interfaces entre les éléments d’où le problème de régularité. Zienkiewicz
et Zhu [Zienkiewicz & Zhu 1987] ont proposé d’interpoler les contraintes avec les mêmes
fonctions de forme que le déplacement afin d’obtenir des contraintes continues sur tout le
domaine. On a alors :
σ ∗ = N · σ̄ ∗ (2.3.33)

avec σ̄ ∗ : contraintes aux noeuds et N : fonctions de forme du déplacement.

La clef de ce type d’estimateur est de construire les contraintes aux noeuds
∗
σ̄ . Zienkiewicz et Zhu ont proposé deux méthodes : la méthode de projec-
tion [Zienkiewicz & Zhu 1987] et la méthode ’Superconvergent Patch Recovery’
[Zienkiewicz & Zhu 1992a, Zienkiewicz & Zhu 1992b].

Méthode de projection :
La méthode de projection proposée dans [Zienkiewicz & Zhu 1987] consiste à construire σ̄ ∗
au sens des moindres carrés à partir de σh .
σ̄ ∗ doit minimiser Z
J(σ̄ ) = (N · σ̄ ∗ − σh )2 dΩ
∗
(2.3.34)
Ω
La minimisation donne : Z
N T · (N · σ̄ ∗ − σh ) dΩ = 0 (2.3.35)
Ω
D’où :
Z
∗ −1
σ̄ = A · N T · σh dΩ (2.3.36)
Ω
Z
avec A = NT · N d Ω
Ω
Cet estimateur d’erreur est simple à mettre en oeuvre et moins coûteux que les autres
estimateurs. Ainsworth et al. [Ainsworth et al. 1989] ont montré que cette méthode est très
efficace. Même quand il existe des singularités, ce type d’estimateur fonctionne bien. Il a
été implémenté dans plusieurs codes industriels comme Abaqus [Dassault_Systèmes 2016] ,
Cast3M [CAST3M 2016] et Code_Aster [Code_Aster 2016].

Méthode Superconvergent Patch Recovery (SPR) :

Plus tard, Zienkiewicz et Zhu [Zienkiewicz & Zhu 1992a, Zienkiewicz & Zhu 1992b] ont amé-
lioré leur technique de calcul des contraintes aux noeuds en utilisant des points dit "supercon-
vergents" (points où la solution éléments finis possède des propriété de superconvergence) et
la notion de patch.
L’approche de projection a été améliorée par l’utilisation des propriétés de superconver-
gence de la solution éléments finis en certains points et par la récupération d’une solution
améliorée dans un patch local.
48 Chapitre 2. Raffinement de maillage en mécanique des solides

Un patch est centré sur un noeud du maillage et contient un ensemble de points super-
convergents (voir figure 2.14). Sur chaque patch ω, la méthode SPR consiste à écrire les
contraintes au noeud comme une expansion polynomiale :

σ̄ ∗ = P · a (2.3.37)

avec P la base de l’expansion polynomiale, a = (a1 , a2 , · · · , aP +1 ) le vecteur contenant les

(P+1) coefficients du polynôme dans le patch ω.
On a, par exemple, P = [1, x, y] pour l’expansion linéaire en 2D et P = [1, x, y, x2 , xy, y 2 ]
pour l’expansion quadratique en 2D.

F I G U R E 2.14 – Patch de la méthode SPR [Delmas 2008]

On cherche à minimiser sur chaque patch (au sens des moindres carrés) :
nX
SP R

J(a) = (σh (xi , yi ) − σ̄ ∗ (xi , yi ))2

i=1
nX
SP R

= (σh (xi , yi ) − P (xi , yi ) · a)2 (2.3.38)

i=1

où (xi , yi ) sont des coordonnées des points superconvergents, et nSP R le nombre de points
superconvergents.
La minimisation de (2.3.38) implique que a satisfasse (équation normale) :
nX
SP R nX
SP R

P T (xi , yi )P (xi , yi )a = P T (xi , yi )σh (xi , yi ) (2.3.39)

i=1 i=1
On a :
a = A−1 b (2.3.40)
où
nX
SP R nX
SP R

A= P T (xi , yi )P (xi , yi ) et b = P T (xi , yi )σh (xi , yi )

i=1 i=1
Connaissant a, on peut exprimer les contraintes au noeud σ̄ ∗ = P · a en utilisant les coordon-
nées du noeud au centre de patch. En utilisant les fonctions de forme du déplacement, on
peut construire un champ de contraintes lissées.
2.3. Estimateurs d’erreur a posteriori 49

Ce type d’estimateur est simple à mettre en œuvre. Il est largement étudié et appliqué
à la méthode des éléments finis. Bouillard et al. [Bouillard et al. 1996] ont utilisé ce type
d’estimateur dans le domaine acoustique. Wriggers et Scherf [Wriggers & Scherf 1998] l’ont
appliqué au problème de contact sans frottement. L’étude dans le cadre de la simulation de
forgeage a été réalisée par Boussetta et al. [Boussetta et al. 2006].

Bilan du chapitre
A partir du problème modèle élastostatique linéaire, nous avons présenté les deux princi-
pales familles de méthodes de raffinement local et les trois principales familles d’estimateurs
d’erreur a posteriori.
Dans notre contexte d’étude, nous souhaitons utiliser une méthode de raffinement
de maillage permettant de considérer le solveur comme une "boîte noire". Les méthodes
de raffinement de type multi-grilles locales sont ainsi choisies. Les travaux précédents
[Barbié 2013, Barbié et al. 2014, Barbié et al. 2015] ont montré les performances de la mé-
thode LDC. Nous continuerons d’utiliser cette méthode dans ce mémoire. Afin de détecter
automatiquement les zones du maillage à raffiner, nous décidons de combiner la méthode
LDC avec un estimateur d’erreur a posteriori. Nous choisissons d’utiliser l’estimateur d’erreur
a posteriori de Zienkiewicz et Zhu (ZZ) [Zienkiewicz & Zhu 1987] en raison de sa simplicité
de mise en oeuvre, de sa rapidité et de son efficacité démontrée dans la littérature.
Chapitre 3

Stratégie de raffinement automatique

de maillage

Après avoir présenté différentes méthodes de raffinement local et différents estimateurs

d’erreur a posteriori dans le Chapitre 2, nous allons proposer ici une stratégie de raffinement
automatique de maillage basée sur la méthode multi-grilles "Local Defect Correction" et l’esti-
mateur a posteriori de Zienkiewicz et Zhu. Nous vérifierons l’efficacité de notre stratégie sur
des cas test étudiant la réponse élastique de la gaine soumise à des chargements simulant le
contact avec la pastille.

Sommaire
3.1 Méthode “Local defect correction” . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.1.1 Formulation de la méthode LDC . . . . . . . . . . . . . . . . . . . . . . . 54
3.1.2 Application la méthode LDC dans le cadre de la méthode des éléments finis 57
3.2 Combinaison de la méthode LDC et de l’estimateur ZZ . . . . . . . . . . . . . 60
3.3 Présentation des cas tests 2D issus de la simulation de l’IPG . . . . . . . . . . 63
3.3.1 Mise en diabolo : cas 2D(r,z) . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3.2 Fissuration de la pastille : cas 2D(r,θ) . . . . . . . . . . . . . . . . . . . . 64
3.4 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.1 Choix numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.2 Erreur relative en norme énergie . . . . . . . . . . . . . . . . . . . . . . 65
3.4.3 Erreur absolue en norme infinie . . . . . . . . . . . . . . . . . . . . . . . 77
3.1. Méthode “Local defect correction” 53

3.1 Méthode “Local defect correction”

Les problèmes industriels peuvent conduire à des solutions possédant des singularité lo-
calisées. La taille de maille dans ces zones nécessite d’être suffisamment petite pour que la
solution atteigne une précision exigée.
Malgré l’intérêt des maillages uniformes (conduisant à des systèmes linéaires pour lesquels
des solveurs rapides et efficaces existent), la taille de maille requise par ces zones de singula-
rité locale rend généralement l’utilisation de tels maillages trop coûteuse.
On utilise alors souvent des maillages non uniformes raffinés localement (voir la section
2.2.1) : dans une ou plusieurs zones, le maillage a un plus grand raffinement afin que la
solution ait une précision équivalente partout. Cependant, la convergence de ces méthodes
est souvent difficile car le système devient mal conditionné.
Une alternative réside dans les méthodes multi-grilles locales qui permettent de conserver les
avantages des maillage uniformes tout en raffinant localement le maillage dans les zones de
singularité.
Nous nous intéressons ici à l’apport de ce type de méthode pour simuler précisément l’inter-
action mécanique pastille-gaine.
Nous nous focaliserons sur la méthode Local Defect Correction (LDC) proposée par Hack-
busch [Hackbusch 1984] car elle se prête bien à une discrétisation par la méthode des élé-
ments finis et a déjà été utilisée avec succès sur des problèmes mécaniques avec singularité lo-
cales dans des travaux précédents [Barbié 2013, Barbié et al. 2014, Barbié et al. 2015]. Nous
en présentons ici les principales caractéristiques.
A partir d’un maillage grossier défini sur l’ensemble du domaine, la méthode LDC consiste à
ajouter récursivement des maillages locaux fins dans des zones où la solution n’est pas suffi-
samment précise (voir figure 3.1). Ces maillages fins peuvent en particulier être uniformes.
La résolution du problème est basée sur un processus itératif qui peut être représenté par des
∧-cycles (voir figure 3.2) : en utilisant des opérateurs de prolongement et de restriction, la
solution du maillage grossier est améliorée par les solutions obtenues sur les maillages fins
locaux. Le processus multi-grilles est obtenu par une généralisation du processus bi-grilles
(deux niveaux de grille) qui peut être résumé ainsi :
• résolution sur le maillage grossier
• définition des conditions aux limites sur le maillage fin à partir de la solution grossière
par un opérateur de prolongement
• résolution sur le maillage fin
• transport des informations du maillage fin vers le maillage grossier par un opérateur
de restriction
• correction de la solution grossière par ajout du résidu calculé à partir de la solution fine
restreinte sur la zone de raffinement.
Ce processus itératif est réalisé jusqu’à la convergence de la solution grossière.
Dans la section suivante, nous détaillons l’algorithme multi-grilles LDC et dérivons son
54 Chapitre 3. Stratégie de raffinement automatique de maillage

Niveau 2
Grille fine Gl∗

Niveau 1

Grille grossière G0

Initialisation Prolongement (conditions aux limites)

Niveau 0

Lissage ou résolution exacte Restriction (correction)

Solution convergée
F I G U R E 3.1 – Exemple de
maillages locaux fins hiérar- F I G U R E 3.2 – Représentation du processus multi-grille lo-
chiques cal pour 4 niveaux de grilles [Barbié 2013]

application à la méthode des éléments finis. Nous nous inspirons des travaux de Hackbush
[Hackbusch 1984], de [Link] et al. [Khadra et al. 1996] et de L. Barbié [Barbié 2013].

3.1.1 Formulation de la méthode LDC

Afin de formuler la méthode LDC, nous considérons un problème défini sur le domaine
borné Ω ⊂ Rn , n ∈ N∗ :
(
Lu = f dans Ω
(P) (3.1.1)
C.L. naturelles sur Γ
avec L un opérateur différentiel en espace et f le second membre. Les "C.L. naturelles" sont
des conditions aux limites choisies pour avoir existence et unicité de la solution.
Nous définissons alors des domaines imbriqués Ωl , 0 ≤ l ≤ l∗ , Ωl ⊂ Ωl−1 avec Ω0 = Ω. Sur
chaque domaine local Ωl , on cherche à résoudre le problème :
(
L|Ωl u = f|Ωl dans Ωl
(P|Ωl ) (3.1.2)
C.L. données par l’algorithme LDC sur Γl

On discrétise chaque domaine Ωl sur grille Gl de taille de maille hl avec hl < hl−1 et de
frontière Γl = ΓD,l ∪ ΓN,l . La solution approchée ul sur Gl est la solution du problème noté
(Pl ) :
(Pl ) : Ll ul = fl + gl dans Gl (3.1.3)

où Ll est la discrétisation de l’opérateur différentiel L|Ωl , fl est la discrétisation du second

membre f|Ωl , gl est la discrétisation des conditions aux limites.
La méthode LDC consiste à résoudre plusieurs fois ce problème (Pl ) avec des seconds
membres corrigés. Ainsi au k ième ∧- cycle, (Pl ) s’écrit :

(Plk ) : Ll ukl = flk + glk dans Gl (3.1.4)

3.1. Méthode “Local defect correction” 55

Notons que par soucis de clarté, les illustrations sont présentées sur des maillages
hiérarchiques avec un rapport de raffinement égal à 2. Cependant, la méthode LDC
est très générique, et les sous-grilles locales ne sont pas forcément hiérarchiques (voir
[Nefedov & Mattheij 2002]).

[Link] Initialisation

Nous notons fl0 la discrétisation sur Gl du second membre f|Ωl du problème continu. La
solution grossière initiale u00 est obtenue en résolvant le problème (P00 ) défini sur tout le
domaine avec une grille G0 .

[Link] Étape de prolongement

Pour résoudre le problème Plk , k ≥ 1, 1 ≤ l ≤ l∗ sur la grille Gl , on doit définir des

conditions aux limites sur la frontière Γl (cf. figure 3.3) :
• S’il existe une frontière commune entre Γl et Γ, on utilise les conditions aux limites
naturelles du problème (P).
• Sur la partie de frontière Γl \Γ 6= ∅, on applique des conditions aux limites de type
Dirichlet (la choix de conditions de type Dirichlet va assurer l’existence et l’unicité de
la solution raffinée). Les valeurs de ces conditions aux limites sont obtenues via un
l
opérateur de prolongement Pl−1 appliqué à la solution du niveau immédiatement plus
k
grossier ul−1 .
ukl|Γl \Γ = (Pl−1
l
ukl−1 )|Γl \Γ (3.1.5)

Conditions aux limites naturelles

Gl
Conditions aux limites de Dirichlet
obenues par prolongation

Gl−1

F I G U R E 3.3 – Prolongement des conditions aux limites

Une fois que les conditions aux limites sur la grille Gl sont définies, la solution ukl est
calculée en résolvant le problème (Plk ).
56 Chapitre 3. Stratégie de raffinement automatique de maillage

[Link] Étape de restriction

La solution ukl , k ≥ 1, l∗ > l ≥ 0 est corrigée via un résidu calculé à partir de la solution
du niveau immédiatement plus fin ukl+1 . Les conditions aux limites définies durant l’étape de
prolongement précédente sont conservées afin de résoudre le nouveau problème (Plk ).
Deux ensembles de noeuds de Gl doivent être définis (voir figure 3.4 et
[Hackbusch 1984]). On introduit Al l’ensemble des noeuds de Gl strictement inclus dans
Gl+1 , et Ål l’ensemble des noeuds Al ne faisant intervenir que des noeuds de Al ∪ (Γl ∩ Γ)
dans le schéma de discrétisation.
Tout d’abord, on restreint la solution ukl+1 sur l’ensemble des noeuds de Al :

ũkl (x) = Rl+1

l
ukl+1 (x) ∀x ∈ Al (3.1.6)
l
où Rl+1 est un opérateur de restriction de Gl+1 vers Gl .
Le résidu sur la grille Gl est alors calculé sur les noeuds de Ål :

rlk (u) = Ll ũkl (x) − fl0 ∀x ∈ Ål (3.1.7)

(a) Zone de restriction Al (b) Zone de correction Ål (exemple pour

un opérateur discrétisé sur un stencil à 5
points)

F I G U R E 3.4 – Zone de restriction et zone de correction

On corrige le second membre :

flk = fl0 + χÅl rlk (u) (3.1.8)

où χÅl est la fonction caractéristique de Ål :

(
1 si x ∈ Ål
χÅl = (3.1.9)
0 si x ∈
/ Ål

On obtient ukl en résolvant le problème (Plk ) avec le second membre corrigé.

Vu que la correction n’est effectuée que sur la zone Ål , Hackbusch [Hackbusch 1984] a
montré que la zone de raffinement doit être suffisamment large pour que la méthode LDC soit
efficace.
3.1. Méthode “Local defect correction” 57

Remarque : L’équation (3.1.8) suppose que la résolution du problème fin soit suffisamment
précise. Dans ce cas, la restriction du résidu du problème fin r̂lk = R̂l+1
l
(Ll+1 (ukl+1 )−fl+1
k k
−gl+1 )
l
(où R̂l+1 est un opérateur de restriction basé sur la conservation des forces) est négligeable
devant rlk :
kr̂lk k ≪ krlk k (3.1.10)

où k · k désigne une norme (choix arbitraire).

[Link] Choix des opérateurs

Le choix des opérateurs de prolongation et de restriction est un point important de l’ef-

ficacité de la méthode LDC. Ils doivent avoir un ordre cohérent avec l’ordre du schéma de
discrétisation afin de conserver l’ordre optimal attendu de convergence. En pratique, on choisit
souvent des interpolateurs d’ordre désiré.

[Link] Algorithme LDC multi-niveaux

Les étapes [Link] à [Link] conduisent à l’algorithme LDC multi-niveaux décrit dans
l’algorithme 1 en page 58.

3.1.2 Application la méthode LDC dans le cadre de la méthode des éléments

finis
Nous appliquons ici la méthode LDC (section [Link]) à une discrétisation par éléments
finis dans le cas particulier de l’élasticité linéaire.
Le problème (Plk ) s’écrit alors sous forme matricielle :

[Kl ][Ulk ] = [Flk ] + [Gl ]Γl ∩ΓN

(
(3.1.11)
[Hl ][Ulk ] = [Ul,D
k
]Γl \Γ + [Ul,D ]Γl ∩ΓD

où [Kl ] est la matrice de rigidité calculée sur la grille Gl , [Flk ] est la matrice des forces
externes volumiques, [Gl ]Γl ∩ΓN est la matrice des forces externes surfaciques (conditions de
Neumann), notée [Gl ] dans la suite, [Hl ] est la matrice de discrétisation des conditions aux
k
limites de Dirichlet, [Ul,D ]Γl \Γ est la matrice des valeurs de Dirichlet interpolées sur la frontière
k
fictive Γl \Γ, notée [Ul,D ] dans la suite, et [Ul,D ]Γl ∩ΓD est la matrice des valeurs de Dirichlet
provenant des conditions aux limites naturelles, notée [Ul,D ] dans la suite.
L’algorithme 2 en page 59 décrit la méthode LDC dans ce cas.
Remarque : Sur le niveau le plus grossier, k=0, seuls les efforts volumiques F0k vont varier
au cours des itérations (invariance des conditions aux limites).
58 Chapitre 3. Stratégie de raffinement automatique de maillage

Algorithme 1 : Algorithme LDC multi-niveaux

Initialisation : calcul de u00 solution du problème initial (P00 );
/*Itération : calcul de ukl */;
pour k=1 à k ∗ faire
uk0 = uk−1
0 ;
/*Étape de prolongement*/;
pour l = 1 à l∗ faire
si k=1 alors
Définition des conditions aux limites naturelles sur Γl ∩ Γ à partir de celles définies
sur Γ ;
Définition du second membre fl0 sur Gl à partir de f
fin
Détermination des conditions aux limites de Dirichlet sur Γl \Γ

ukΓl \Γ = (Pl−1
l
ukl−1 )|Γl \Γ

Calcul de ukl en résolvant du problème (Plk ) avec flk = flk−1 ;

fin
/*Étape de restriction*/;
pour l = l∗ − 1 à 0 faire
Restriction de la solution fine ukl+1 sur Al
l
ũkl (x) = Rl+1 ukl+1 (x) ∀x ∈ Al

Calcul du résidu local grossier sur Ål

rlk (u) = Ll ũkl (x) − fl0 ∀x ∈ Ål

Correction du second membre

flk = fl0 + χÅl rlk (u)

Calcul de ukl en résolvant du problème (Plk )

fin
fin
3.1. Méthode “Local defect correction” 59

Algorithme 2 : Algorithme LDC appliquée à la méthode des éléments finis

Initialisation : calcul de la solution [U00 ] ;
(
[K0 ][U00 ] = [F00 ] + [G0 ]
[H0 ][U00 ] = [U0,D ] sur ΓD,0

/*Itération : calcul de [Ulk ]*/;

pour k=1 à k ∗ faire
[U0k ] = [U0k−1 ];
/*Étape de prolongement*/;
pour l = 1 à l∗ faire
si k=1 alors
Définition des conditions aux limites de Dirichlet [Ul,D ] et de Neumann [Gl ] sur
Γl ∩ Γ à partir de celles définies sur Γ;
Définition du seconde membre [Fl0 ] par discrétisation f sur Gl ;
fin
Détermination des conditions aux limites sur Γl \Γ;

k l k
[Ul,D ] = [Pl−1 ][Ul−1 ] sur Γl \Γ
Calcul de [Ulk ];
(
[Kl ][Ulk ] = [Flk−1 ] + [Gl ]
[Hl ][Ulk ] = [Ul,D
k
] + [Ul,D ] sur ΓD,l

fin
/*Restriction : Correction sur les grilles grossières*/;
pour l = l∗ − 1 à 0 faire
k
Restriction de la solution fine [Ul+1 ] sur Al

[Ũlk ] = [Rl+1
l k
][Ul+1 ] sur Al

Calcul du résidu local sur Ål

[Dlk ] = [Kl ][Ũlk ] − [Fl0 ] sur Ål

Correction du second membre

[Flk ] = [Fl0 ] + χÅl [Dlk ]

Calcul de [Ulk ] en résolvant :

(
[Kl ][Ulk ] = [Flk ] + [Gl ]
[Hl ][Ulk ] = [Ul,D
k
] + [Ul,D ] sur ΓD,l

fin
fin
60 Chapitre 3. Stratégie de raffinement automatique de maillage

3.2 Combinaison de la méthode LDC et de l’estimateur ZZ

Afin de détecter automatiquement les zones du maillage à raffiner, nous décidons de
combiner la méthode LDC avec un estimateur d’erreur a posteriori. Comme présenté dans la
section 2.3, nous choisissons d’utiliser l’estimateur d’erreur a posteriori de Zienkiewicz et Zhu
(ZZ) [Zienkiewicz & Zhu 1987] en raison de sa simplicité de mise en oeuvre, de sa rapidité et
de son efficacité démontrée dans la littérature. L’idée ici est de générer automatiquement des
sous-niveaux indépendemment de la connaissance de l’utilisateur et du problème traité.
La méthode LDC étant très générique, elle se combine facilement avec un estimateur d’er-
reur a posteriori. La combinaison de la méthode LDC et de l’estimateur ZZ est réalisée de façon
séquentielle. Comme dans les travaux de Barbié et al. [Barbié 2013, Barbié et al. 2014], nous
appliquons uniquement l’estimateur d’erreur durant la première étape de prolongement afin
de générer les différents niveaux de grilles fines. Pendant la première étape de prolongement,
après chaque résolution, l’estimateur ZZ est appliqué sur ce niveau afin de détecter la zone à
raffiner qui va générer une nouvelle grille. Les grilles ainsi générées seront conservées pen-
dant tout le processus LDC. Ce processus génère donc automatiquement les grilles imbriquées
Gl , 1 ≤ l ≤ l ∗ .
Dans la stratégie proposée par Barbié et al. [Barbié et al. 2014], pour chaque grille Gl ,
les auteurs ont choisi de raffiner les éléments dont l’indicateur donné par l’estimateur ZZ
en norme énergie relative était supérieur à α fois la valeur maximale de cet indicateur avec
α une constante donnée entre 0 et 1.
Cette stratégie est simple à implémenter et s’applique à n’importe quel estimateur d’erreur
et n’importe quelle norme. Les résultats numériques dans [Barbié et al. 2014] ont montré
l’efficacité de cette méthode.
Cependant, les zones détectées sont dépendantes du choix du coefficient α. Ce coefficient est
difficile à trouver a priori et dépendant du problème. Par définition, il y a toujours des éléments
qui vérifient la condition de raffinement. Un critère d’arrêt doit donc être ajouté afin de
stopper la génération de sous-niveaux dans la méthode LDC. Barbié et al. [Barbié et al. 2014]
ont proposé que ce critère d’arrêt s’appuie sur un nombre minimal d’éléments dans la zone à
raffiner.
En vue d’améliorer la façon d’automatiser la détection des zones à raffiner et ainsi la
procédure de raffinement, nous proposons ici de raffiner tous les éléments possédant une
erreur estimée en contraintes supérieure d’un seuil fixé par l’utilisateur. C’est-à-dire que nous
raffinons les éléments K de la grille Gl pour lesquels l’erreur estimée locale eK vérifie :
keK k > β (3.2.1)
où k · k désigne une norme choisie par l’utilisateur et β le seuil d’erreur demandé par l’utilisa-
teur.
Remarque : Dans le cas où l’on veut conserver des maillages réglés, il faudra ajouter
des éléments supplémentaires afin d’obtenir des maillages locaux réglés. Par ailleurs, cette
structuration de maillage rend l’estimateur ZZ plus fiable (cf. coin rentrant). Afin de garantir
que tous les noeuds détectés soient dans l’ensemble A0l (voir section [Link]), nous élargissons
la zone à raffiner d’une rangée d’éléments.
3.2. Combinaison de la méthode LDC et de l’estimateur ZZ 61

L’algorithme 3 en page 62 décrit le couplage de la methode LDC avec l’estimateur d’erreur

a posteriori de Zienkiewicz et Zhu.

Critère de convergence du processus itératif :

Le processus itératif (∧-cycle) est considéré comme convergé si l’erreur relative entre les deux
solutions grossières successives en norme L2 est plus petite que δ :

kuk0 − uk−1
0 kL2 (Ω)
≤δ (3.2.2)
kuk0 kL2 (Ω)

La méthode LDC converge souvent très vite. Typiquement, pour δ = 10−5 , moins de quatre
itérations sont nécessaires pour que la solution converge.

Critère supplémentaire d’arrêt de génération des sous-niveaux

La condition (3.2.1) sert également de condition d’arrêt. Quand les indicateurs d’erreur de
tous les éléments sont inférieurs au seuil β, le processus de génération de sous-niveaux s’ar-
rête. Cependant, quand le problème présente des singularités, ce critère est difficile voire
impossible à atteindre et le raffinement continu indéfiniment. Nous proposons ici des critères
supplémentaires d’arrêt pour la génération des grilles locales.
Dans un premier temps, nous utiliserons le critère présenté dans les travaux de Barbié
et al. [Barbié et al. 2014]. Il consiste à définir un nombre d’éléments minimal conduisant à
générer un sous-niveau supplémentaire. Cependant, ce nombre d’éléments est arbitrairement
choisi alors qu’il semble devoir dépendre du problème traité.
Nous proposons de comparer ce critère à un critère géométrique. Le critère géométrique
s’appuie sur la mesure de la zone détectée par l’équation (3.2.1) et est lié aux limites de la
méthode LDC : quand la taille du maillage le plus fin est trop petite, la solution grossière n’est
plus impactée par la solution la plus fine. Nous proposons deux types de critère géométrique.
Le premier critère consiste à arrêter le raffinement lorsque la mesure de la zone détectée est
inférieure à la mesure d’un élément de maillage grossier de départ (Sl < Se voir figure 3.5).
S’il existe des parties connexes détectées, chacune est traitée de manière indépendante vis-à-
vis du critère d’arrêt. Le second critère est basé sur le ratio de surface de la zone à raffiner
vis-à-vis du domaine initial. Si ce ratio est inférieur à ǫ (Sl < ǫS0 voir figure 3.5), la génération
des sous-niveaux s’arrête.

Sl

S0

Se

F I G U R E 3.5 – Critère d’arrêt géométrique · S0 mesure de surface de maillage grossier de

départ · Sl mesure de surface de la zone détectée · Se mesure d’un élément du maillage
grossier de départ
62 Chapitre 3. Stratégie de raffinement automatique de maillage

Algorithme 3 : Algorithme LDC appliquée à la méthode des éléments finis couplée avec
l’estimateur d’erreur
Initialisation : Calcul de la solution [U00 ] ;
(
[K0 ][U00 ] = [F00 ] + [G0 ]
[H0 ][U00 ] = [U0,D ] sur ΓD,0

/*Itération : calcul de [Ulk ]*/;

pour k=1 à k ∗ faire
[U0k ] = [U0k−1 ];
/*Étape de prolongement*/;
pour l = 1 à l∗ faire
si k=1 alors
/*Génération des sous-niveaux et définition des conditions aux limites sur Γl ∩ Γ*/;
1
Application de l’estimateur d’erreur sur la solution [Ul−1 ] afin de détecter la zone à
raffiner;
Raffinement de la zone détectée : obtention de la grille Gl ;
Définition des conditions aux limites de Dirichlet [Ul,D ] et de Neumann [Gl ] sur
Γl ∩ Γ à partir de celles définies sur Γ;
Définition du seconde membre [Fl0 ] par discrétisation f sur Gl ;
fin
Détermination des conditions aux limites sur Γl \Γ
k l k
[Ul,D ] = [Pl−1 ][Ul−1 ] sur Γl \Γ

Calcul de [Ulk ] (
[Kl ][Ulk ] = [Flk−1 ] + [Gl ]
[Hl ][Ulk ] = [Ul,D
k
] + [Ul,D ] sur ΓD,l

fin
/*Étape de restriction*/;
pour l = l∗ − 1 à 0 faire
k
Restriction de la solution fine [Ul+1 ] sur Al

[Ũlk ] = [Rl+1
l k
][Ul+1 ] sur Al

Calcul du résidu local sur Ål

[Dlk ] = [Kl ][Ũlk ] − [Fl0 ] sur Ål

Correction du second membre

[Flk ] = [Fl0 ] + χÅl [Dlk ]

Calcul de [Ulk ] en résolvant :

(
[Kl ][Ulk ] = [Flk ] + [Gl ]
[Hl ][Ulk ] = [Ul,D
k
] + [Ul,D ] sur ΓD,l

fin
fin
3.3. Présentation des cas tests 2D issus de la simulation de l’IPG 63

3.3 Présentation des cas tests 2D issus de la simulation de l’IPG

Nous avons proposé une stratégie de raffinement automatique de maillage combinant
la méthode LDC et l’estimateur ZZ avec plusieurs critères d’arrêt dans section 3.2. Nous
présentons dans cette section des cas tests 2D simplifiés issus de la simulations de l’IPG pour
évaluer les performances de cette stratégie.
Les cas tests utilisés ici sont inspirés de Marelle et Thouvenin [Marelle & Thouvenin 2011].
Nous n’étudions que la réponse de la gaine ayant un comportement élastique avec un module
d’Young E = 100 GPa et un coefficient de Poisson ν = 0, 3. Le contact entre la pastille
et la gaine est représenté par une pression discontinue sur le rayon interne de la gaine.
Nous présentons ici deux cas tests (2D axisymétrie et 2D déformation plane) afin d’étudier
séparément les deux phénomènes caractéristiques de l’Interaction mécanique Pastille-Gaine
(IPG) : la déformation en diabolo et la fissuration de la pastille (voir la section 1.4).

3.3.1 Mise en diabolo : cas 2D(r,z)

Le cas test 2D axisymétrique, appelé par la suite cas test 2D(r,z), représente l’effet du
phénomène de mise en diabolo de la pastille sur la gaine (voir section 1.3.1 et section 1.4). Le
contact avec la pastille est modélisé par un pic de pression sur 600 µm au dessus du plan inter-
pastille. Pour des raisons de symétrie, seule une demi-hauteur de pastille est modélisée (voir
figure 3.6). Sur la figure 3.6b, l’ensemble des conditions aux limites permettant de résoudre
le problème est donné.
4.1mm 0.6mm

Conditions de symétrie
(Relation d’ensemble)
0.6 mm
Plan Médian-Pastille

Pression interne constante

due à la pastille (80MPa)
13.5 mm

Pression externe
6.15 mm

constante (15.5MPa)

Pic de pression
0.6 mm

due à la forme de
diabolo (150MPa) Plan Inter-Pastille
Conditions de symétrie

(a) Déformation en diabolo (b) Définition du problème

F I G U R E 3.6 – Cas test 2D(r,z)

Ce cas test a une géométrie très simple : des maillages réguliers structurés uniformes
peuvent parfaitement représenter le domaine d’étude.
Comme il n’existe pas de solution analytique pour ce cas test, nous utilisons une solution
64 Chapitre 3. Stratégie de raffinement automatique de maillage

obtenue avec un maillage très fin (maillage uniforme avec un pas de maille d’environ 2 µm et
1053801 nœuds en total ) comme solution de référence.

3.3.2 Fissuration de la pastille : cas 2D(r,θ)

Le second cas test, appelé cas test 2D(r,θ) par la suite, respecte l’hypothèse des déforma-
tions planes. Ce cas test est destiné à simuler l’effet de la fissuration de la pastille sur la gaine.
La pastille est supposée se fissurer de façon régulière (voir [Nonon et al. 2004] et figure 3.7a).
Pour des raisons de symétrie, on modélise seulement 1/16 de la gaine (voir figure 3.7b). L’ou-
verture de la fissure est de 8 µm. La fissure de la pastille est représentée par une discontinuité
de pression sur la gaine. La définition du problème traité est donnée sur la figure (3.7b).

métrie
de sy
itions o
Cond 22,5

Pression externe
Pression interne constante (15.5MPa)
constante (80MPa)

Discontinuité de
pression (0MPa) 8 µm
R = 4,2 mm 600 µm
Conditions de symétrie

(a) Fissuration de pastille (b) Définition du problème

F I G U R E 3.7 – Cas test 2D(r,θ)

Pour ce cas test, une solution analytique a été proposée par Roberts [Roberts 1978]. Ce-
pendant, cette solution est exprimée sous la forme d’une série Fourier et le nombre de termes
requis afin d’obtenir une solution assez précise devient très vite prohibitif. Comme pour le
cas test 2D(r,z), on utilise une solution obtenue avec un maillage très fin (maillage quasi-
uniforme avec un pas de maille d’environ 1 µm et 260365 nœuds en total) comme solution
de référence.

3.4 Résultats numériques

Nous présentons dans cette section les résultats numériques de notre stratégie de raffine-
ment (voir la section 3.2) avec les différents critères d’arrêt de génération des sous-niveaux
pour différentes normes : erreur relative en norme énergie, erreur absolue en norme infinie.

3.4.1 Choix numériques

Seules des opérations de pré et post traitement sont effectuées pour appliquer l’algorithme
LDC combiné avec l’estimateur ZZ. Le solveur éléments finis peut être vu comme une boîte
noire donc la stratégie proposée ici est applicable à n’importe quel solveur éléments finis. Le
3.4. Résultats numériques 65

logiciel Cast3M [CAST3M 2016] a été utilisé ici. Nous choisissons d’utiliser des éléments Q1
de Lagrange standard car le logiciel Cast3M ne traite pas correctement le problème de contact
avec des éléments supérieurs à 1 même pour les éléments Pk .
Les éléments finis Q1 permettent théoriquement d’atteindre une convergence d’ordre 2 en
norme L2 (cf. section 2.1.4). Cependant, l’approximation à l’ordre 1 de la discontinuité de
pression par le maillage conduit notre schéma de discrétisation à être du premier ordre
[Ramière 2008, Barbié et al. 2014]. Des opérateurs de prolongement et de restriction d’ordre
1 sont suffisants. Nous choisissons d’utiliser comme opérateur de prolongement l’interpolateur
bilinéaire et comme opérateur de restriction une restriction canonique (car on utilise des
maillage hiérarchiques).
À chaque étape, nous résolvons les problèmes par un solveur direct. À la première prolonga-
tion, la factorisation de la matrice de rigidité est conservée pour les résolutions des ∧-cycles
suivants. Pour les problèmes linéaires, seuls les temps de résolutions de la première étape de
prolongment sont donc coûteux.
Les calculs sont effectués sur un ordinateur avec 8 processeurs de Intel Core i7-4790 et un
espace de mémoire de 8 Gigabytes.

3.4.2 Erreur relative en norme énergie

[Link] Généralités

Nous commençons par étudier les résultats obtenus lorsque l’erreur choisie dans le cri-
tère (3.2.1) est l’erreur relative et la norme est la norme énergie.
L’erreur relative en norme énergie pour l’élément K donnée par l’estimateur ZZ peut être
exprimée sous la forme suivante :

∗ 1/2
− σh ) : (ε∗ − εh ) dK
R
K (σ
keK kE = R
∗ ∗
(3.4.1)
K σ : ε dK

où σh représente les contraintes calculées par la méthode des éléments finis, σ ∗ les contraintes
lissées, C le tenseur d’élasticité et ε = C −1 σ, ε∗ = C −1 σ ∗ .
On peut démontrer que lorsque la génération des sous-niveaux s’arrête grâce au critère
keK kE ≤ β avec β seuil d’erreur alors l’erreur keΩ kE est également inférieure au seuil, avec
keΩ kE l’erreur sur le maillage composite (réunion des sous-niveaux) .

∗ 1/2
− σh ) : (ε∗ − εh ) dK
R
K (σ
keK kE = R
∗ ∗
≤β
K σ : ε dK
∗ σ ) : (ε∗ − εh ) dK
R
K (σ −R h
⇐⇒ ∗ : ε∗ dK
≤ β2
K σ
Z Z
⇐⇒ ∗ ∗ 2
(σ − σh ) : (ε − εh ) dK ≤ β ( σ ∗ : ε∗ dK)
K K
66 Chapitre 3. Stratégie de raffinement automatique de maillage

En sommant sur le maillage composite, on a :

∗ +1
lX Z ∗ +1
lX Z
X X
∗ ∗ 2
=⇒ (σ − σh ) : (ε − εh ) dK ≤ β ( σ ∗ : ε∗ dK)
l=1 K⊂Ωl−1 \Ωl K l=1 K⊂Ωl−1 \Ωl K

avec Ωl∗ +1 = ∅.

On a donc :
 ∗ +1
1/2
lX X Z
 (σ ∗ − σh ) : (ε∗ − εh ) dK 
K⊂Ωl−1 \Ωl K
 
 l=1 
keΩ kE =  ∗
 ≤β (3.4.2)
 lX +1 X Z 
∗ ∗
 
 σ : ε dK 
l=1 K⊂Ωl−1 \Ωl K

Ainsi, les erreurs composites entre la solution de la méthode LDC et la solution de référence
à la fin du calcul (après convergence des ∧-cycles) nous permettrons de vérifier si l’estimateur
ZZ est un bon estimateur (i.e. s’il estime bien keK kE ).

[Link] Résultats pour le cas test 2D(r,z)

Nous appliquons la stratégie de couplage LDC-ZZ présentée dans la section 3.2 au cas test
2D(r,z).
On utilise un rapport de raffinement égal à 2. Les sous-grilles générées à partir d’un
maillage initial de pas 164 µm, pour un seuil d’erreur β fixé à 1% et un critère d’arrêt de 10
éléments minimaux sont présentées sur la figure 3.8. A partir du calcul sur le maillage initial,
on estime les erreurs relatives sur chaque élément avec l’estimateur ZZ. Nous choisissons
alors les éléments possédant une erreur supérieure au seuil de 1% (éléments en vert) pour
générer un sous-niveau. Nous ajoutons ainsi récursivement des sous-niveaux jusqu’à atteindre
le critère d’arrêt. Ces premiers résultats semblent très satisfaisants puisque les sous-grilles se
localisent automatiquement autour de la zone de discontinuité de pression. On peut en effet
penser que c’est dans cette zone qu’il est nécessaire d’avoir un maillage assez fin pour avoir
une solution précise.

Critère d’arrêt en nombre d’élément

Le couplage LDC-ZZ a été testé pour différents pas de maillages de départ (328 µm,164 µm,
82 µm, 41 µm) et différents seuils d’erreur β (5%, 2%, 1%). Les résultats sont présentés
dans le tableau 3.1 pour le critère d’arrêt en nombre d’éléments (NM IN = 10). Le résultat
pour chaque cas test est présenté en quatre lignes : la première ligne indique le nombre
de sous-niveaux générés automatiquement, la deuxième ligne contient le nombre de nœuds
de chaque niveau, la troisième ligne présente l’erreur relative entre la solution composite à
convergence et la solution de référence et la quatrième ligne monte le temps total de calcul.
3.4. Résultats numériques 67

SCAL SCAL SCAL SCAL

< 3.55E−02 < 2.44E−02 < 1.86E−02 < 1.60E−02
> 6.88E−03 > 2.55E−03 > 8.93E−04 > 1.05E−03

3.00E−02 3.00E−02 3.00E−02 3.00E−02

2.50E−02 2.50E−02 2.50E−02 2.50E−02

2.00E−02 2.00E−02 2.00E−02 2.00E−02

1.50E−02 1.50E−02 1.50E−02 1.50E−02

1.00E−02 1.00E−02 1.00E−02 1.00E−02

8.20E−03 8.20E−03 8.20E−03 8.20E−03

6.40E−03 6.40E−03 6.40E−03 6.40E−03

4.60E−03 4.60E−03 4.60E−03 4.60E−03

2.80E−03 2.80E−03 2.80E−03 2.80E−03

1.00E−03 1.00E−03 1.00E−03 1.00E−03

Maillage grossier 1er sous-niveau 2e sous-niveau 3e sous-niveau

F I G U R E 3.8 – Génération des sous-niveaux · Cas test 2D(r,z) · Maillage initial de pas 164 µm
et seuil d’erreur relative en norme énergie des contraintes β est fixé à 1% · Critère d’arrêt : 10
éléments. Pour chaque niveau de maillage, les erreurs élémentaires estimées sont présentées
à gauche et les mailles en vert présentés à droite sont les éléments à raffiner.

Nous remarquons que pour certains seuils, les maillages sont complètement raffinés. Dans ce
cas, les maillages totalement raffinés deviennent alors le nouveau maillage de départ.
La première conclusion qui peut être tirée du tableau 3.1 est que la génération des sous-
niveaux s’arrête automatiquement pour le nombre minimal d’éléments choisi. Quelque soit le
maillage de départ, les zones raffinées sont les mêmes (taille et nombre d’éléments) pour un
seuil fixé. Il montre aussi l’efficacité de la méthode LDC. En effet, quelque soit le maillage de
départ, on obtient le même niveau d’erreur pour un pas de sous-niveau le plus fin équivalent.
Les erreurs diminuent rapidement en ajoutant des sous-niveaux locaux avec un nombre limité
de degrés de liberté. Il semble donc plus intéressant en terme de nombre de nœuds de com-
mencer avec un maillage grossier et d’ajouter des sous-niveaux locaux que de commencer avec
un maillage très fin. La seconde conclusion que l’on peut tirer est que la stratégie proposée
fonctionne bien : à la fin des calculs, les erreurs relatives obtenues sont inférieures aux seuils
demandés. De plus, chaque sous-niveau généré (excepté pour le maillage hi/8) est généra-
lement nécessaire pour que l’erreur reste inférieure au seuil demandé. Ce cas test présente
une géométrie simple avec petit nombre de degrés de liberté sur laquelle des maillages réglés
uniformes peuvent être appliqués. Les temps d’exécution restent très faibles même quand on
rajoute des sous-niveaux. Les temps de calculs pour les trois premiers maillages initiaux sont
comparables. Par contre, les temps pour le maillage de départ de hi/8 sont plus élevés. La
tendance qui se dégage en terme de temps de calcul est qu’il y a donc un compromis à trouver
68 Chapitre 3. Stratégie de raffinement automatique de maillage

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 1 2 4
nombre de noeuds par niveau : 125/99 125/315/697 125/441/1105/2343/247
hi
erreur relative : 4,84% 2,11% 0,653%
temps de calcul : 0.1 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 1 3
nombre de noeuds par niveau : 441 441/697 441/1105/2343/247
hi/2
erreur relative : 4,20% 2,01% 0,653%
temps de calcul : 0.03 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 0 2
nombre de noeuds par niveau : 1649 1649 1649/2343/247
hi/4
erreur relative : 1,92% 1,92% 0,550%
temps de calcul : 0.2 s 0.2 s 1.3 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 6369 6369 6369/247
hi/8
erreur relative : 0,757% 0,757% 0,443%
temps de calcul : 1.7 s 1.7 s 3.0 s

Tableau 3.1 – Cas test 2D(r,z) pour différents pas de maillages de départ et différents seuils
fixés avec la formule (3.4.1) · hi=328 µm · Critère d’arrêt : 10 éléments

entre le maillage de départ et le nombre de sous-niveaux.

Pour le cas test avec le pas de maillage de départ le plus fin (hi/8) et pour le seuil de 1%,
le nombre de sous-niveaux générés est supérieur au nombre minimal de niveaux requis pour
respecter ce seuil. Il y a deux raisons possibles à cela :
— Soit les erreurs données par ZZ surestiment l’erreur de référence.
— Soit la condition keK k ≤ β est une condition suffisante mais pas nécessaire.
Pour vérifier la première hypothèse, nous avons comparé les erreurs estimées par l’estima-
teur de ZZ et les erreurs obtenues sur le même maillage entre la solution éléments finis et la
solution de référence. La figure 3.9 montre que les erreurs sont en réalité sous-estimées par
l’estimateur d’erreur ZZ.
Il semblerait donc que la condition keK kE ≤ β soit suffisante mais non nécessaire pour
que keΩ kE ≤ β. Pour assurer que l’erreur globale soit inférieure au seuil, on n’a pas donc
besoin que l’erreur de chaque élément soit inférieure au seuil. En effet, l’erreur relative globale
est une sorte de moyenne des erreurs relatives locales.
Une autre version de l’estimateur ZZ en norme énergie est introduite par
[Zienkiewicz & Zhu 1987] afin d’assurer que le dénominateur de la formule (3.4.1) soit
non nul : R ∗ ∗ 1/2
K (σ −Rσh ) : (ε − εh ) dK

keK kE = R
∗ ∗
(3.4.3)
K σh : εh dK + K (σ − σh ) : (ε − εh ) dK
Cette formule d’erreur (3.4.3) est la formule utilisée dans les codes industriels comme
3.4. Résultats numériques 69

SCAL SCAL
< 3.54E−02 < 9.18E−02
> 6.88E−03 > 7.22E−03

9.00E−02 9.00E−02
8.00E−02 8.00E−02
7.00E−02 7.00E−02
6.00E−02 6.00E−02
5.00E−02 5.00E−02
4.00E−02 4.00E−02
3.00E−02 3.00E−02
2.00E−02 2.00E−02
1.00E−02 1.00E−02
8.20E−03 8.20E−03
6.40E−03 6.40E−03
4.60E−03 4.60E−03
2.80E−03 2.80E−03
1.00E−03 1.00E−03

Erreurs par rapport à

Erreurs estimées ZZ la solution référence

F I G U R E 3.9 – Comparaison entre les erreurs estimées et les erreurs calculées par rapport à la
solution de référence · Cas test 2D(r,z) · Pas de maillage 164 µm.
70 Chapitre 3. Stratégie de raffinement automatique de maillage

Cast3M, Code_Aster, Catia, Abaqus etc. Nous testons les mêmes seuils et les mêmes maillages
de départ avec cette formule d’erreur (3.4.3). Les résultats sont présentés dans le tableau 3.2.

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 1 2 4
nombre de noeuds par niveau : 125/99 125/315/697 125/441/1105/2343/247
hi
erreur relative : 4,84% 2,11% 0,653%
temps de calcul : 0.1 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 1 3
nombre de noeuds par niveau : 441 441/697 441/1105/2343/247
hi/2
erreur relative : 4,20% 2.01% 0,653%
temps de calcul : 0.03 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 0 2
nombre de noeuds par niveau : 1649 1649 1649/2343/247
hi/4
erreur relative : 1,92% 1,92% 0,550%
temps de calcul : 0.2 s 0.2 s 1.3 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 6369 6369 6369/247
hi/8
erreur relative : 0,757% 0,757% 0,443%
temps de calcul : 1.7s s 1.7 s 3.0 s

Tableau 3.2 – Cas test 2D(r,z) pour différents pas de maillages de départ et seuils fixés avec la
formule (3.4.3) · hi=328 µm · Critère d’arrêt : 10 éléments

On remarque que les résultats sont similaires à ceux obtenus avec la formule (3.4.1). Les
deux formules sont donc équivalentes pour nos cas test. Dans la suite, la formule (3.4.3) sera
utilisée pour rester compatible avec codes industriels.
Les tableaux 3.1 et 3.2 montrent que l’application au cas test 2D(r,z) de la stratégie de
raffinement LDC-ZZ est prometteuse. Cependant, nous avons utilisé le critère d’arrêt basé
sur un nombre d’éléments arbitraire qu’il est difficile de déterminer à priori. La génération
des sous-niveaux peut ne pas s’arrêter pour un autre nombre d’éléments choisi et ce nombre
d’éléments varie en fonction du cas test.

Critères d’arrêt géométriques

Nous allons tester maintenant les critères d’arrêt géométriques proposés dans la section 3.2.
Nous commençons par le critère d’arrêt basé sur la mesure de la zone détectée qui doit
être inférieure à la mesure d’un élément de maillage grossier de départ. Les résultats sont
présentés dans le tableau 3.3.
Les résultats montrent que ce critère fonctionne correctement. L’arrêt de la génération des
sous-niveaux semble plus optimal, notamment pour les maillages initiaux les plus grossiers.
Le principal inconvénient de cette stratégie est qu’elle dépend de la taille de maille initiale.
3.4. Résultats numériques 71

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 1 2 3
nombre de noeuds par niveau : 125/99 125/315/697 125/441/1105/1155
hi
erreur relative : 4,84% 2,11% 0,969%
temps de calcul : 0.1 s 0.3 s 0.7 s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 441 441/697 441/1105/2343
hi/2
erreur relative : 4,20% 2.01% 0,899%
temps de calcul : 0.03 s 0.2 s 0.8s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 1649 1649/117 1649/2343/247
hi/4
erreur relative : 1,92% 1,01% 0,550%
temps de calcul : 0.2 s 0.3 s 1.3 s
nombre de sous-niveaux : 0 0 2
nombre de noeuds par niveau : 6369 6369 6369/247/165
hi/8
erreur relative : 0,757% 0,757% 0,440%
temps de calcul : 1.7 s 1.7 s 3.2 s

Tableau 3.3 – Cas test 2D(r,z) pour différents pas de maillages de départ et seuils fixés avec la
formule (3.4.3) · hi=328 µm · Critère d’arrêt : mesure inférieure à un élément du maillage
de départ

Quand la taille de maille initiale est petite, cette stratégie conduit à générer des sous-niveaux
non nécessaires (cf. pas de maillage de départ hi/8).
Nous allons tester le critère d’arrêt basé sur le ratio de surface de la zone à raffiner du
domaine de départ. Ici, nous choisissons 0,5% de la mesure du domaine de départ comme
critère d’arrêt. Les résultats sont présentés dans le tableau 3.4. Ces résultats sont satisfaisants.
En effet, ce critère permet d’arrêter la génération des sous-niveaux indépendamment de la
taille de maille initiale. Le nombre de sous-niveaux ajoutés est suffisant pour assurer que
les erreurs relatives obtenues sont inférieures aux seuils demandés. Cependant, il conduit à
générer un sous-niveau de plus que nécessaire pour les calculs avec le seuil β = 1%

Génération de sous-niveaux avec un algorithme de type FMG local

Nous avons essayé d’utiliser l’algorithme FMG local (voir section 2.2.2 et figure 3.10) afin
de générer un sous-niveau supplémentaire uniquement lorsque l’erreur composite globale
estimée après convergence était supérieure au seuil demandé.
Les résultats pour différents pas de maillages de départ et différents seuils sont présentés
dans le tableau 3.5. Le résultat pour chaque cas test est présenté en cinq lignes : la première
ligne indique le nombre de sous-niveaux générés automatiquement, la deuxième ligne contient
le nombre de nœuds de chaque niveau, la troisième ligne donne l’erreur relative estimée sur
la solution composite à convergence, la quatrième ligne présente l’erreur relative entre la
72 Chapitre 3. Stratégie de raffinement automatique de maillage

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 1 3 4
nombre de noeuds par niveau : 125/99 125/315/697/117 125/441/1105/2343/247
hi
erreur relative : 4,84% 1,34% 0,653%
temps de calcul : 0.1 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 2 3
nombre de noeuds par niveau : 441 441/697/117 441/1105/2343/247
hi/2
erreur relative : 4,20% 1,17% 0,653%
temps de calcul : 0.04 s 0.3 s 1.2 s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 1649 1649/117 1649/2343/247
hi/4
erreur relative : 1,92% 1,02% 0,550%
temps de calcul : 0.1 s 0.2 s 1.3 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 6369 6369 6369/247
hi/8
erreur relative : 0,757% 0,757% 0,443%
temps de calcul : 1.7 s 1.7 s 3.0 s

Tableau 3.4 – Cas test 2D(r,z) pour différents pas de maillages de départ et seuils fixés avec
la formule (3.4.3) · hi=328 µm · Critère d’arrêt : mesure inférieure à 0,5% de la mesure du
domaine initial

Grille fine Gl∗

Grille grossière G0

Lissage Restriction du résidu

Résolution exacte Prolongement de la correction

Solution convergée

F I G U R E 3.10 – Représentation du cycle FMG pour 4 niveaux de grilles

solution composite à convergence et la solution de référence et la cinquième ligne montre le

temps de calcul. Comme les erreurs sont sous-estimées par l’estimateur ZZ (voir figure 3.12),
cette stratégie conduit à ne pas suffisamment raffiner. Les erreurs par rapport à la solution de
référence dépassent le seuil à la fin du calcul. L’algorithme FMG combiné à l’estimateur ZZ ne
permet pas d’obtenir une stratégie de raffinement plus optimale que les ∧-cycles standards de
la méthode LDC.
Par ailleurs, du fait que les sous-niveaux de la méthode LDC sont générés à la première
3.4. Résultats numériques 73

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 125 125/315 125/441/1105
hi erreur relative estimée : 3,83% 1,99% 0,986%
erreur relative : 8,52% 4,25% 1,95%
temps de calcul : 0.01 s 0.2 s 0.4 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 441 441 441/1105
hi/2 erreur relative estimée : 1,89% 1,89% 0,986%
erreur relative : 4,20% 4,20% 1,95%
temps de calcul : 0.03 s 0.03 s 0.2 s
nombre de sous-niveaux : 0 0 0
nombre de noeuds par niveau : 1649 1649 1649
hi/4 erreur relative estimée : 0,924% 0.924% 0,924%
erreur relative : 1,92% 1,92% 1,92%
temps de calcul : 0.2 s 0.2 s 0.2 s
nombre de sous-niveaux : 0 0 0
nombre de noeuds par niveau : 6369 6369 6369
hi/8 erreur relative estimée : 0,455% 0,455% 0,445%
erreur relative : 0,757% 0,757% 0,757%
temps de calcul : 1.7 s 1.7 s 1.7 s

Tableau 3.5 – Cas test 2D(r,z) pour différents pas de maillages de départ et seuils fixés avec
l’algorithme FMG local · hi =328 µm · Critère d’arrêt : erreur composite estimée inférieure au
seuil

étape de prolongement et donc à partir de solutions non convergées, nous nous demandons si
une procédure de détection du deuxième niveau après convergence sur le premier sous-niveau
(convergence du ∧-cycle bi-grille) peut améliorer la stratégie proposée.
En comparant les erreurs estimées sur le premier sous-niveau après convergence du cycle
bi-grille (cf. figure 3.11a) et celles estimées sur la solution non convergée à la première étape
de prolongement (cf. figure 3.11b), nous pouvons conclure que les zones détectées restent
toujours les mêmes. Le processus FMG ne permet donc pas vraiment d’améliorer la stratégie
LDC-ZZ proposée.

Conclusion
Nous avons appliqué plusieurs stratégies de génération de sous-niveaux au cas test 2D(r,z).
La stratégie de raffinement LDC-ZZ et les critères d’arrêt proposés permettent tous de générer
des sous-niveaux et de stopper la génération de façon automatique. À la fin des calculs, les
erreurs par rapport à la solution de référence sont généralement inférieures aux seuils. La
stratégie basée sur l’algorithme FMG local ne permet pas d’améliorer les solutions obtenues
74 Chapitre 3. Stratégie de raffinement automatique de maillage

SCAL SCAL
> 2.56E−03 > 2.56E−03
< 2.45E−02 < 2.44E−02

2.43E−02 2.42E−02
2.32E−02 2.32E−02
2.22E−02 2.22E−02
2.12E−02 2.11E−02
2.01E−02 2.01E−02
1.91E−02 1.90E−02
1.80E−02 1.80E−02
1.70E−02 1.70E−02
1.59E−02 1.59E−02
1.49E−02 1.49E−02
1.39E−02 1.38E−02
1.28E−02 1.28E−02
1.18E−02 1.18E−02
1.07E−02 1.07E−02
9.69E−03 9.67E−03
8.65E−03 8.63E−03
7.60E−03 7.59E−03
6.56E−03 6.55E−03
5.52E−03 5.51E−03
4.48E−03 4.47E−03
3.43E−03 3.43E−03

(a) après convergence (b) premier prolongement

F I G U R E 3.11 – Erreurs estimées du premier sous-niveau après convergence à gauche et du

premier prolongement à droite · Cas test 2D(r,z) · Maillage initial de pas 164 µm · Seuil de
contraintes fixé à 1%.

car l’erreur est trop sous-estimée par l’estimateur ZZ.

Parmi les différents critères d’arrêt, les résultats obtenus avec le critère d’arrêt en nombre
d’élément sont satisfaisants. Les erreurs par rapport à la solution de référence sont en général
inférieures aux seuils. Cependant, le nombre d’éléments minimal est difficile à définir et varie
en fonction du problème. Les résultats avec le critère d’arrêt lié à la mesure d’un élément de
maillage grossier de départ sont optimaux pour les maillages de départ de hi et hi/2 avec
le seuil β = 1%. Cependant, pour les maillages de départ de hi/4 et hi/8, des sous-niveaux
inutiles sont générés. Les résultats pour le critère d’arrêt basé sur le ratio de surface de la zone
à raffiner vis-à-vis du domaine de départ montrent que les erreurs obtenues sont inférieures
aux seuils. Il semble tous les sous-niveaux pour les seuils des 5% et 2% sont nécessaires.
Cependant, quelques sous-niveaux inutiles pour le seuil de 1% sont ajoutés pour les maillages
les plus grossiers.

[Link] Résultats du cas test 2D(r,θ)

Pour aider à faire un choix entre les critères d’arrêt indépendant du problème, nous testons
les différents critères d’arrêt de génération des sous-niveaux sur le cas test 2D(r,θ). Les résul-
3.4. Résultats numériques 75

tats sont présentés dans le tableau 3.6 pour le critère en nombre d’élément (ici, le nombre
optimal semble être 7), dans le tableau 3.7 pour le critère géométrique de la mesure d’un
élément du maillage de départ et dans le tableau 3.8 pour le critère géométrique de du ratio
de la mesure du domaine initial.

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 0 0 2
nombre de noeuds par niveau : 231 231 231/861/121
hi
erreur relative : 7,05% 7,05% 1,27%
temps de calcul : 0.01 s 0.01 s 0.3 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 861 861 861/121
hi/2
erreur relative : 3,19% 3,19% 1,27%
temps de calcul : 0.08 s 0.08 s 0.1 s
nombre de sous-niveaux : 0 0 0
nombre de noeuds par niveau : 3321 3321 3321
hi/4
erreur relative : 1,22% 1,22% 1,22%
temps de calcul : 0.5 s 0.5 s 0.5 s

Tableau 3.6 – Cas test 2D(r,θ) pour différents pas de maillages de départ et seuils fixés avec
la formule (3.4.3) · hi=109 µm · Critère d’arrêt : 7 éléments

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 231 231/81 231/861/121
hi
erreur relative : 7,05% 3,26 % 1,27%
temps de calcul : 0.03 s 0.1 s 0.3 s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 861 861/81 861/121/121
hi/2
erreur relative : 3,19 % 1,27% 0,462%
temps de calcul : 0.07 s 0.2 s 0.3 s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 3321 3321/81 3321/121/121
hi/4
erreur relative : 1,22% 0,294% 0,289%
temps de calcul : 0.8 s 1.0 s 1.2 s

Tableau 3.7 – Cas test 2D(r,θ) pour différents pas de maillages de départ et seuils fixés avec
la formule (3.4.3) · hi=109 µm · Critère d’arrêt : mesure inférieure à un élément de maillage
de départ
76 Chapitre 3. Stratégie de raffinement automatique de maillage

Seuil d’erreur β
Maillage de départ 5% 2% 1%
nombre de sous-niveaux : 0 2 3
nombre de noeuds par niveau : 231 231/99/81 231/861/121/121
hi
erreur relative : 7,05% 1,48 % 0,464%
temps de calcul : 0.02 s 0.02 s 0.4 s
nombre de sous-niveaux : 0 1 2
nombre de noeuds par niveau : 861 861/81 861/121/121
hi/2
erreur relative : 3,19 % 1,27% 0,462%
temps de calcul : 0.07 s 0.08 s 0.5 s
nombre de sous-niveaux : 0 0 1
nombre de noeuds par niveau : 3321 3321 3321/121
hi/4
erreur relative : 1,22% 1,22% 0,296%
temps de calcul : 0.6 s 0.6 s 1.0 s

Tableau 3.8 – Cas test 2D(r,θ) pour différents pas de maillages de départ et seuils fixés avec
la formule (3.4.3) · hi=109 µm · Critère d’arrêt : mesure inférieure 0,5% de la mesure du
domaine initial

Les résultats montrent que la génération des sous-niveaux s’arrête bien automatiquement
dans tous les cas. La combinaison LDC-ZZ fonctionne bien ici aussi puisqu’à la fin des calculs,
les erreurs obtenues sont généralement inférieures aux seuils fixés. Le nombre d’éléments
minimal pour arrêter la génération de sous-niveaux a été difficile à définir pour le critère d’ar-
rêt en nombre d’élément. En effet, le cas test 2D(r,θ) est plus sensible au nombre d’éléments
minimal que le cas test 2D(r,z). Pour le maillage de départ de hi, il a 3 éléments possédant une
erreur supérieure à 2%. Afin qu’un sous-niveaux soit ajouté, le nombre minimal d’élément doit
être inférieur à 3. Cependant, pour le seuil à 1%, la génération de sous-niveaux ne s’arrête
pas si le nombre minimal d’élément est inférieur à 7. Le deuxième critère améliore la stratégie
pour les maillages de départ avec hi et hi/2 mais on voit bien ici que ça conduit à raffiner
inutilement pour hi/4. Ce critère d’arrêt n’est pas non plus optimal. Les résultats 3.8 montrent
que le critère géométrique basé sur le rapport de mesure par rapport au domaine initial est le
meilleur pour ce cas test. Chaque sous-niveau généré semble être nécessaire pour que l’erreur
de la solution LDC par rapport à la solution de référence reste inférieure au seuil.
Les résultats obtenus pour le maillage le plus grossier (avec le pas de maillage hi = 109µm)
ne sont pas aussi bons qu’attendu pour les deux premiers seuils (5% et 2%), surtout pour
les deux premiers critères d’arrêt. En vue d’en trouver la raison, l’erreur relative estimée par
l’estimateur de Zienkiewicz et Zhu et l’erreur relative par rapport à la solution de référence
sont comparées dans la figure 3.12. La figure 3.12 montre que l’erreur la plus importante
par rapport à la solution de référence est située autour du coin en bas à droite du maillage.
Ceci est dû à la mauvaise approximation des conditions aux limites du fait du maillage trop
grossier. Ce type d’erreur ne peut pas être détecté par l’estimateur ZZ qui est basé sur le lissage
de la solution. Si la solution est lisse mais fausse, l’estimateur ZZ ne peut pas la détecter. Dans
3.4. Résultats numériques 77

cette zone, l’estimateur ZZ sous-estime énormément l’erreur.

SCAL SCAL
> 8.53E−03 > 1.05E−02
< 2.95E−02 < 1.58E−01

0.15 0.15

0.13 0.13

0.11 0.11

9.00E−02 9.00E−02

7.00E−02 7.00E−02

5.00E−02 5.00E−02

3.00E−02 3.00E−02

2.60E−02 2.60E−02

2.20E−02 2.20E−02

1.80E−02 1.80E−02

1.40E−02 1.40E−02

1.00E−02 1.00E−02

F I G U R E 3.12 – Erreur estimée par l’estimateur ZZ (à gauche) et erreur entre la solution EF

et la solution de référence (à droite) (pas de maillage hi = 109 µm)

[Link] Bilan

Nous avons appliqué la stratégie de raffinement et les critères d’arrêt proposés pour les
deux cas tests 2D(r,z) et 2D(r, θ). Nous avons comparé la stratégie proposée de couplage
LDC-ZZ et la stratégie basée sur l’algorithme FMG local. En raison de la sous-estimation de
l’estimateur ZZ, la stratégie de couplage FMG local-ZZ fonctionne moins bien que la stratégie
de couplage LDC-ZZ. Les résultats de la stratégie de couplage LDC-ZZ sont satisfaisants : les
sous-niveaux sont ajoutés automatiquement dans les zones d’intérêt, les erreurs de solutions
diminuent rapidement par ajout de sous-niveaux, les erreurs sont de même niveau pour
des maillages les plus fins équivalents. Parmi les trois critères d’arrêt permettant stopper
automatiquement la génération de sous-niveaux, le critère d’arrêt basé sur le ratio de surface
de la zone à raffiner vis-à-vis du domaine de départ semble être le plus optimal. Dans la suite,
nous allons tester cette stratégie pour un critère de détection basé sur l’erreur absolue en
norme infinie.

3.4.3 Erreur absolue en norme infinie

[Link] Généralités

Dans le cadre de la modélisation de l’IPG, l’erreur absolue en norme infinie sur les
contraintes est souvent beaucoup plus intéressante d’un point de vue pratique pour les méca-
niciens. L’erreur absolue en norme infinie estimée par ZZ peut être exprimée sous la forme
suivante :

keK,abs k∞ = kσ ∗ − σh k∞ (3.4.4)
78 Chapitre 3. Stratégie de raffinement automatique de maillage

Cependant, cette norme est rarement utilisée pour les estimateurs d’erreur a posteriori
dans la littérature car elle repose directement sur la "fiabilité" quantitative de l’estimateur.

[Link] Résultats pour le cas test 2D(r,z)

Dans un premier temps, nous comparons les erreurs absolues estimées par l’estimateur de
ZZ et les erreurs absolues par rapport à la solution de référence pour le cas test 2D(r,z) sur les
contraintes équivalentes de von Mises et sur différentes composantes du champ de contraintes
dans la figure 3.13. Cette figure montre que les erreurs sur la contrainte équivalente de von
Mises et les erreurs sur les contraintes les contraintes cisaillement sont celles qui sont les
mieux estimées par l’estimateur ZZ. De plus, la contrainte équivalente de von Mises prend
en compte toutes les des composantes des contraintes. C’est donc plus intéressant d’estimer
l’erreur absolue en norme infinie de la contrainte équivalente de von Mises.
Nous allons donc nous baser sur les erreurs absolues sur les contraintes équivalentes de
von Mises pour détecter les éléments à raffiner. A la fin du calcul, afin de valider la stratégie
proposée nous évaluerons les erreurs absolues entre la solution de LDC et la solution de
référence sur le maillage composite privé du dernier niveau. En effet, du fait de la singularité
de contrainte, le dernier niveau contient toujours une zone sur laquelle l’erreur maximale est
supérieure au seuil. Les résultats sont présentés dans le tableau 3.9.
Seuil d’erreur β
Maillage de départ 10MPa 5MPa 2MPa
nombre de sous-niveaux : 4 5 6
hi nombre de noeuds par niveau : 125/297/697/891/285 125/441/1071/2541/2925/621 125/441/1649/4653/11895/12513/38293
erreur absolue : 8,34MPa 5,00MPa 2,25MPa
nombre de sous-niveaux : 3 4 5
hi/2 nombre de noeuds par niveau : 441/697/891/285 441/1071/2541/2925/621 441/1649/4653/11895/12513/38239
erreur absolue : 8,32MPa 5,00MPa 2,25MPa
nombre de sous-niveaux : 2 3 4
hi/4 nombre de noeuds par niveau : 1649/891/315 1649/2541/2925/621 1649/4653/11895/12513/38239
erreur absolue : 8,28MPa 4,99MPa 2,25MPa
nombre de sous-niveaux : 1 2 3
hi/8 nombre de noeuds par niveau : 6369/315 6369/2925/621 6369/11895/12513/38239
erreur absolue : 8,29MPa 4,99MPa 2,25MPa

Tableau 3.9 – Cas test 2D(r,z) pour différents pas de maillages de départ et seuils fixés · Erreur
absolue sur la contrainte équivalente de von Mises en norme infinie · hi =328 µm · Critère
d’arrêt : mesure inférieure à 0,5% de la mesure du domaine initial

Les résultats obtenus en utilisant l’erreur absolue sur la contrainte équivalente de von Mises
en norme infinie sont très satisfaisants pour le cas test 2D(r,z). La stratégie LDC-ZZ marche
très bien dans ce cas. Elle permet d’arrêter automatiquement de génération de sous-niveaux.
Les sous-niveaux sont ajoutés dans les même zones et avec le même nombre d’éléments
pour un seuil fixé. Pour la taille de maille la plus fine équivalente, les erreurs absolues sont
similaires. A la fin des calculs, les erreurs absolues en norme infinie entre la solution éléments
finis et la solution de référence sont généralement très proches des erreurs demandées.
3.4. Résultats numériques 79

SCAL SCAL
< 6.02E+07 < 9.99E+07 SMRR SMRR SMZZ SMZZ
>−1.22E+07 >−1.29E+06 < 2.05E+07 < 5.76E+07 < 3.94E+07 < 9.28E+07
>−2.05E+07 >−3.41E+07 >−3.94E+07 >−6.15E+07

9.00E+06 9.00E+06
9.00E+06 9.00E+06 9.00E+06 9.00E+06
8.20E+06 8.20E+06
8.20E+06 8.20E+06 8.20E+06 8.20E+06
7.40E+06 7.40E+06
7.40E+06 7.40E+06 7.40E+06 7.40E+06
6.60E+06 6.60E+06
6.60E+06 6.60E+06 6.60E+06 6.60E+06
5.80E+06 5.80E+06
5.80E+06 5.80E+06 5.80E+06 5.80E+06
5.00E+06 5.00E+06
5.00E+06 5.00E+06 5.00E+06 5.00E+06
4.20E+06 4.20E+06
4.20E+06 4.20E+06 4.20E+06 4.20E+06
3.40E+06 3.40E+06
3.40E+06 3.40E+06 3.40E+06 3.40E+06
2.60E+06 2.60E+06
2.60E+06 2.60E+06 2.60E+06 2.60E+06
1.80E+06 1.80E+06
1.80E+06 1.80E+06 1.80E+06 1.80E+06
1.00E+06 1.00E+06
1.00E+06 1.00E+06 1.00E+06 1.00E+06
ZZ-VonMises Erreur référence
-VonMises
ZZ-SMRR Erreur référence-SMRR ZZ-SMRZ Erreur référence-SMRZ
(a) Contrainte équiva-
(b) σRR (c) σZZ
lente de von Mises
SMTT SMTT SMRZ SMRZ
< 1.69E+07 < 7.03E+06 < 3.30E+07 < 2.87E+07
>−1.69E+07 >−5.83E+07 >−3.30E+07 >−5.12E+07

9.00E+06 9.00E+06 9.00E+06 9.00E+06

8.20E+06 8.20E+06 8.20E+06 8.20E+06

7.40E+06 7.40E+06 7.40E+06 7.40E+06

6.60E+06 6.60E+06 6.60E+06 6.60E+06

5.80E+06 5.80E+06 5.80E+06 5.80E+06

5.00E+06 5.00E+06 5.00E+06 5.00E+06

4.20E+06 4.20E+06 4.20E+06 4.20E+06

3.40E+06 3.40E+06 3.40E+06 3.40E+06

2.60E+06 2.60E+06 2.60E+06 2.60E+06

1.80E+06 1.80E+06 1.80E+06 1.80E+06

1.00E+06 1.00E+06 1.00E+06 1.00E+06

ZZ−SMTT Erreur référence-SMTT ZZ-SMRZ Erreur référence-SMRZ

(d) σT T (e) σRZ

F I G U R E 3.13 – Comparaison des erreurs absolues estimées par l’estimateur ZZ et des erreurs
absolues entre la solution des éléments finis et la solution de référence · Cas test 2D(r,z) hi
=164 µm

[Link] Résultats pour le cas test 2D(r,θ)

Afin de conforter les conclusions précédentes, nous appliquons la stratégie de couplage

LDC-ZZ basée sur les erreurs absolues en norme infinie sur les contraintes équivalentes de von
Mises pour le cas test 2D(r,θ) avec différents pas de maillages de départ et différents seuils
80 Chapitre 3. Stratégie de raffinement automatique de maillage

(10 MPa, 5 MPa,2 MPa). Les résultat sont présentés dans le tableau 3.10.

Seuil d’erreur β
Maillage de départ 10MPa 5MPa 2MPa
nombre de sous-niveaux : 3 3 3
hi nombre de noeuds par niveau : 231/99/99/99 231/861/143/169 231/861/3321/289
erreur absolue : 4,45MPa 3,35MPa 2,06MPa
nombre de sous-niveaux : 2 2 2
hi/2 nombre de noeuds par niveau : 861/99/99 861/143/169 861/3321/289
erreur absolue : 2,97MPa 2,80MPa 2,06MPa
nombre de sous-niveaux : 1 1 1
hi/4 nombre de noeuds par niveau : 3321/99 3321/143 3321/289
erreur absolue : 3,04MPa 2,62MPa 2,06MPa

Tableau 3.10 – Cas test 2D(r,θ) pour différents pas de maillages de départ et seuils fixés ·
Erreur absolue sur la contrainte équivalente de von Mises en norme infinie · hi=109 µm ·
Critère d’arrêt : mesure inférieure à 0,5% de la mesure du domaine initial

Les résultats montrent que la stratégie de couplage LDC-ZZ utilisant l’erreur absolue sur
la contrainte équivalente de von Mises en norme infinie fonctionne également bien pour
le cas test 2D(r,θ). La génération des sous-niveaux s’arrête automatiquement et de manière
cohérente. Les sous-niveaux sont ajoutés dans les mêmes zones et ont les mêmes nombres
d’éléments pour des calculs avec le même seuil et différents pas de maillages de départ. Les
erreurs en norme absolue diminuent rapidement par ajout des sous-niveaux. Les erreurs par
rapport à la solution de référence après la convergence de la méthode LDC sont proches pour
des tailles de maille fines équivalentes. Cependant, il semble que cette stratégie génère plus de
sous-niveaux que nécessaire pour les cas tests avec seuils fixé à 10 MPa. Nous allons regarder
en détail l’évolution de l’erreur en fonction de nombre de sous-niveaux. Pour ce faire, nous
reprenons les cas test avec le seuil fixé à 10 MPa et différents pas de maillages de départ. Pour
chaque cas test, nous fixons le nombre de sous-niveaux dans la phase de génération des sous-
niveaux. Les erreurs calculées entre les solutions convergées et la solution de référence sont
présentées dans le tableau 3.11 ci-après. On voit que l’ajout d’un sous-niveau fait descendre
d’un facteur supérieur à 2 l’erreur absolue. Ce tableau montre que les sous-niveaux générés
par la stratégie de raffinement automatique sont en fait nécessaires pour que l’erreur absolue
reste inférieure au seuil.
1 Sous-niveau 2 Sous-niveaux 3 Sous-niveaux
hi 25,73MPa 10,63MPa 4,45MPa
hi/2 10,44MPa 2,97MPa
hi/4 3,04MPa

Tableau 3.11 – Erreur des sous-niveaux pour différents pas de maillages de départ et seuil fixé
à 10MPa · hi=109 µm

Nous pouvons donc conclure que la stratégie de couplage LDC-ZZ avec critère d’arrêt
3.4. Résultats numériques 81

automatique fonctionne globalement très bien pour le cas test 2D(r,θ) en utilisant l’erreur
absolue sur la contrainte équivalente de von Mises en norme infinie.

Bilan du chapitre
Dans le cadre de la mise en place d’une méthode multi-grilles permettant de simuler
l’Interaction mécanique Pastille-Gaine (IPG) de façon précise, nous avons étudié dans ce
chapitre les performances d’une stratégie basée sur le couplage de la méthode multi-grilles
LDC et de l’estimateur d’erreur a posteriori ZZ.
Afin de générer les sous-niveaux automatiquement et indépendamment du problème traité,
nous avons proposé une stratégie de couplage LDC-ZZ basée directement sur l’erreur estimée
par l’estimateur ZZ. Divers critères d’arrêt supplémentaires ont également été proposés afin
d’arrêter la génération des sous-niveaux quand le problème présente des singularités de
contraintes. Les résultats numériques des cas tests 2D élastiques nous ont permis d’apprécier
les avantages et les limitations de chacun des critères d’arrêt. De manière générale, le couplage
LDC-ZZ fonctionne très bien puisque les sous-niveaux sont ajoutés automatiquement autour
des zones de singularités. La stratégie proposée conduit généralement à des erreurs (en norme
énergie ou infinie) inférieures au seuil demandé mais un des critères d’arrêt est plus optimal
en terme de sous-niveaux et d’indépendance au problème traité. Il s’agit du critère d’arrêt
basé sur le ratio de surface de la zone à raffiner sur le domaine initial le plus optimal pour
des cas tests 2D élastiques.
Cette stratégie de couplage LDC-ZZ et ce critère d’arrêt ont été appliqués et ont fonctionné
correctement sur différents cas test 2D élastiques avec différentes normes d’erreur. Cette
stratégie de raffinement semble donc très générique.
Chapitre 4

Application de la stratégie de
raffinement automatique au problème
de contact

L’objectif du chapitre est de tester la stratégie de raffinement automatique de maillage

introduite au chapitre 3 sur un problème de contact. L’enjeu ici est de vérifier que la stratégie
proposée permet un raffinement multi-corps avec contact efficace.

Sommaire
4.1 Problème de contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.1.1 Problème de contact unilatéral . . . . . . . . . . . . . . . . . . . . . . . 86
4.1.2 Problème de contact avec frottement . . . . . . . . . . . . . . . . . . . . 88
4.2 Quelques méthodes pour résoudre le problème de contact avec ou sans
frottement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.2.1 Méthode de pénalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.2.2 Méthode des multiplicateurs de Lagrange . . . . . . . . . . . . . . . . . 92
4.2.3 Méthode du Lagrangien augmenté . . . . . . . . . . . . . . . . . . . . . 92
4.2.4 Méthode de point fixe . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.2.5 Algorithmes de résolutions . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.3 Stratégie de raffinement automatique du maillage adaptée au problème de
contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.4 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.4.1 Présentation du cas test . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.4.2 Illustration de la détection des zones à raffiner . . . . . . . . . . . . . . 97
4.4.3 Résultats numérique pour un cas test sans jeu ni frottement . . . . . . . 98
4.4.4 Résultats numériques pour un cas test avec jeu sans frottement . . . . . 105
4.4.5 Résultats numériques pour un cas test avec jeu avec frottement . . . . . 106
4.4.6 Contact avec frontière courbe . . . . . . . . . . . . . . . . . . . . . . . . 110
4.1. Problème de contact 85

4.1 Problème de contact

Nous considérons un problème statique de contact entre deux solides déformables notés
Ω1 et Ω2 , Ω = Ω1 ∪ Ω2 (voir figure 4.1). Les matériaux de ces deux solides sont supposés être
élastiques linéaires.

F I G U R E 4.1 – Problème modèle F I G U R E 4.2 – Zone de contact

de contact entre deux solides défor-
mables

Les frontières suffisamment régulières de ces deux solides, notées respectivement ∂Ω1 et
∂Ω2 , sont divisées en trois parties ΓD , ΓN , ΓC et ∂Ω = ∂Ω1 ∪∂Ω2 = ΓD ∪ΓN ∪ΓC , ΓD ∩ΓN = ∅ :
— la partie ΓD est la zone où les déplacements sont imposés (conditions aux limites de
Dirichlet).
— la partie ΓN est la zone où les forces sont appliquées (conditions aux limites de Neu-
mann).
— la partie ΓC est la zone où les deux solides sont supposés être en contact. Les deux
parties Γ1C et Γ2C sont les zones de contact respectivement de Ω1 et Ω2 . On a ΓC =Γ1C ∪
Γ2C . A priori, la partie ΓC est inconnue.
La projection de la force surfacique sur les bords de contact ΓC sur la direction normale
et la direction tangentielle donne :

FCN = FC · n
(4.1.1)
FCT = FC − FCN n

avec n le vecteur normal unitaire dirigé vers l’extérieur du domaine.

Remarque : FCN est un scalaire et FCT est un vecteur.
Nous nous intéressons ici aux problèmes de contact unilatéral avec et sans frottement.
Pour plus de détails, le lecteur pourra se reporter à [Duvaut & Lions 1976, Wriggers 2006].
86 Chapitre 4. Application de la stratégie au problème de contact

4.1.1 Problème de contact unilatéral

[Link] Formulation du problème

Nous supposons qu’il n’y a ni pénétration ni adhésion entre les deux solides dans la zone
de contact. Le contact unilatéral sans frottement se traduit par une composante tangentielle
nulle des forces de contact. Les forces de contact sont donc normales.
Sur la zone de contact, la condition de non-pénétration revient à imposer que le déplace-
ment relatif uN entre Γ1C et Γ2C est inférieur au jeu initial d :

u N = u 1 · n1 + u 2 · n2 ≤ d (4.1.2)

avec u1 , u2 les déplacements respectivement sur Γ1C , Γ2C et n1 , n2 sont les vecteurs unitaires
normaux dirigés vers l’extérieur de Γ1C , Γ2C (voir figure 4.2). On a ainsi n1 = −n2 sous
l’hypothèse des petites perturbations.
La condition de non-adhésion impose que les efforts normaux ne peuvent correspondre
qu’à de la compression, c’est-à-dire la force normale de contact est toujours non-positive :

FCN ≤ 0 (4.1.3)

Il faut noter qu’il n’existe alors que deux statuts entre les deux solides : soit contact
(uN = d, FCN ≤ 0), soit décollement (uN ≤ d, FCN = 0). Nous obtenons alors la condition de
complémentarité :
(uN − d)FCN = 0 (4.1.4)

La loi de contact unilatéral (4.1.2), (4.1.3) et (4.1.4) peut être présentée par la fi-
gure 4.3. Sous l’hypothèse des petites perturbations, le problème élastostatique avec ces
conditions conduit au problème de Signorini (voir [Duvaut & Lions 1976, Raous et al. 1988,
Lebon & Raous 1992]) qui peut s’exprimer sous la forme :


 −div σ = f dans Ω



 σ = Cε dans Ω
1


T



 ε= grad u + grad u dans Ω
 2
u = u0 sur ΓD (4.1.5)





 σn = FN sur ΓN
N N
uN ≤ d, FC ≤ 0, (uN − d)FC = 0




 sur ΓC
N T T
σn = FC n + FC , FC = 0 sur ΓC


[Link] Formulation variationnelle

Nous définissons l’espace des champs cinématiquement admissibles :

V = {v ∈ (H 1 (Ω))m | v = u0 sur ΓD } (4.1.6)

avec m dimension de l’espace.

4.1. Problème de contact 87

décollement
jeu

contact

F I G U R E 4.3 – Loi de contact

Soit VC le sous-ensemble de V contenant des champs vérifiant les conditions de non-

pénétration (4.1.2) :

1 2
VC = {v ∈ V | vN = vN · n1 + v N · n2 ≤ d sur ΓC } (4.1.7)

VC est un convexe fermé.

La formulation variationnelle en déplacements s’écrit :
Z Z Z
σ : ε(v − u) dΩ = f · (v − u) dΩ + (σn) · (v − u), ∀v ∈ Vc (4.1.8)
Ω Ω ∂Ω

En décomposant les éléments liés à la frontière, on a :

Z Z Z
(Cε(u)) : ε(v − u) dΩ = f · (v − u) dΩ + FN · (v − u) dΓ
Ω Ω ΓN
Z Z (4.1.9)
+ (σn) · (v − u0 ) dΓ + (σn) · (v − u) dΓ
ΓD ΓC
Z
Le terme (σn) · (v − u0 ) dΓ et les forces tangentielles sur les bords ΓC sont nuls. On
ΓD
a: Z Z Z
(Cε(u)) : ε(v − u) dΩ = f · (v − u) dΩ + FN · (v − u) dΓ
Ω Ω ΓN
Z (4.1.10)
+ FCN n · (v − u) dΓ
ΓC

Pour simplifier l’expression et être plus général, on utilise les notations suivantes :
Z
a(u, v) := (Cε(u)) : ε(v) dΩ
ZΩ Z
l(v) := f · v dΩ + FN · v dΓ
Ω ΓN
88 Chapitre 4. Application de la stratégie au problème de contact

avec a un opérateur bilinéaire symétrique, l un opérateur linéaire. On a :

Z
a(u, v − u) − l(v − u) = FCN n · (v − u) dΓ
ZΓC
= FCN (vN − uN ) dΓ (4.1.11)
ΓC
Z
= (FCN (vN − d) − FCN (uN − d)) dΓ
ΓC

Les conditions du problème de Signorini (4.1.4), (4.1.3) et dans l’espace VC nous donnent :

FCN (vN − d) − FCN (uN − d) ≥ 0 (4.1.12)

On obtient donc une inéquation variationnelle :

a(u, v − u) − l(v − u) ≥ 0 (4.1.13)

Du fait de l’inéquation (4.1.13), le problème de contact est non-linéaire (il entre dans le cadre
de l’analyse non régulière).
La formulation variationnelle du problème de Signorini (4.1.5) peut s’écrire sous la forme :
(
Trouver u ∈ VC tel que
(4.1.14)
a(u, v − u) ≥ l(v − u) ∀v ∈ VC

Remarque : on peut démontrer que le problème (4.1.5) est équivalent au problème

(4.1.14).
Duvaut et Lions [Duvaut & Lions 1976] ont montré que la solution du problème (4.1.14)
peut être obtenue par une minimisation de l’énergie potentielle :
(
Trouver u ∈ VC tel que
(4.1.15)
J(u) ≤ J(v) ∀v ∈ VC

avec
1
J(v) = a(v, v) − l(v)
2
Les problèmes (4.1.15) et (4.1.14) sont équivalents. L’existence et d’unicité du problème
de Signorini a été démontrée par Fichera [Fichera 1973].

4.1.2 Problème de contact avec frottement

[Link] Formulation du problème

Nous cherchons une formulation d’un problème de contact avec frottement. Nous utilisons
la loi classique de Coulomb [Duvaut & Lions 1976] qui peut se représenter par les figures 4.4
et 4.5. Il n’existe que deux statuts de frottement :
• statut d’adhérence : si la force tangentielle |FCT | est inférieure à µ|FCN |, alors la vitesse
tangentielle est nulle, u̇T = 0.
4.1. Problème de contact 89

• statut de glissement : si la force tangentielle |FCT | est égale à µ|FCN |, alors la vitesse
tangentielle est non nulle. Elle est colinéaire et de direction opposée à la force de
contact : u̇T = −λFCT .
La loi de Coulomb peut alors s’exprimer sous la forme quasi-statique :

|FCT | ≤ µ|FCN | et
( T
|FC | < µ|FCN | ⇒ u̇T = 0 (adhérence) sur ΓC (4.1.16)
|FCT | = µ|FCN | ⇒ ∃λ > 0, u̇T = −λFCT (glissement)

où | · | désigne la valeur absolue si elle est appliquée à un scalaire et la norme euclidienne

si elle est appliquée à un vecteur, µ est le coefficient de frottement, FCT la force tangentielle
dans la zone de contact ΓC , u̇T la vitesse tangentielle dans la zone de contact ΓC .

glissement adhérence

adhérence glissement glissement

glissement

F I G U R E 4.4 – Loi de Coulomb F I G U R E 4.5 – Cône de Coulomb en 2D

Dans [Raous et al. 1988, Cocu et al. 1996, Raous 1999] est introduite une formulation
"statique" (en déplacement) de la loi de frottement de Coulomb obtenue à partir d’une dis-
crétisation temporelle de la loi en vitesse. À chaque pas de temps, on a alors à résoudre un
problème statique de type (4.1.17). Dans ce mémoire, nous ne considérons que la formulation
statique. Le formulation statique est aussi une approche valide pour les problèmes sous un
chargement monotone.

|FCT | ≤ µ|FCN | et
( T
|FC | < µ|FCN | ⇒ uT = 0 (adhérence) sur ΓC (4.1.17)
|FCT | = µ|FCN | ⇒ ∃λ > 0, uT = −λFCT (glissement)

où uT est le déplacement tangentiel dans la zone de contact ΓC .

90 Chapitre 4. Application de la stratégie au problème de contact

Le problème de Signorini avec frottement de Coulomb peut alors s’exprimer sous la forme :


 −div σ = f dans Ω



 σ = Cε dans Ω
1


grad u + gradT u


ε= dans Ω





 2
u = u0 sur ΓD



σn = FN sur ΓN (4.1.18)

σn = FCN n + FCT sur ΓC




uN ≤ d, FCN ≤ 0, (uN − d)FCN = 0

sur ΓC



 ( T
|FC | < µ|FCN | ⇒ uT = 0



|FCT | ≤ µ|FCN | et


 sur ΓC
|FCT | = µ|FCN | ⇒ ∃λ > 0, uT = −λFCT


[Link] Formulation variationnelle

Nous cherchons la formulation variationnelle du problème de contact avec frottement

(4.1.18). Nous utilisons les mêmes définitions que pour le problème de contact unilatéral
introduites dans la section [Link]. La différence entre les deux cas réside dans le fait que la
force tangentielle de contact n’est pas nulle pour le problème de contact avec frottement. On
a ainsi :
Z Z
N
a(u, v − u) − l(v − u) = FC (vN − uN ) dΓ + FCT · (vT − uT ) dΓ ∀v ∈ VC (4.1.19)
ΓC ΓC

Nous introduisons un terme non différentiable :

Z
j(u, v) = µ|FCN ||vT |dΓ (4.1.20)
ΓC

En utilisant l’inéquation (4.1.12), on a :

Z
a(u, v − u) − l(v − u) + j(u, v) − j(u, u) ≥ FCT · (vT − uT ) dΓ
ΓC
Z Z
N
+ µ|FC ||vT |dΓ − µ|FCN ||uT |dΓ
ΓC ΓC
(4.1.21)
Selon les conditions (4.1.17),
• si |FCT | < µ|FCN |, alors uT = 0, on a :
Z Z Z
T N
FC · (vT − uT ) dΓ + µ|FC ||vT |dΓ − µ|FCN ||uT |dΓ
ΓC ΓC ΓC
Z Z
= FCT · vT dΓ + µ|FCN ||vT |dΓ (4.1.22)
ΓC ΓC
Z
≥ (µ|FCN | − |FCT |)|vT |dΓ ≥ 0
ΓC
4.2. Méthodes pour résoudre le problème de contact 91

• si |FCT | = µ|FCN |, alors uT = −λFCT , on a :

Z Z Z
T N
FC · (vT − uT ) dΓ + µ|FC ||vT |dΓ − µ|FCN ||uT |dΓ
ΓC ΓC ΓC
Z Z Z
= FCT · (vT + λFCT ) dΓ + µ|FCN ||vT |dΓ − µ|FCN ||λFCT |dΓ
Z ΓC Z ΓC ΓC
(4.1.23)
= FCT · vT dΓ + µ|FCN ||vT |dΓ
ΓC ΓC
Z
≥ (µ|FCN | − |FCT |)|vT |dΓ = 0
ΓC

On obtient alors :

a(u, v − u) − l(v − u) + j(u, v) − j(u, u) ≥ 0 (4.1.24)

La formulation variationnelle du problème de Signorini avec frottement de Coulomb

(4.1.18) peut s’écrire sous la forme :
(
Trouver u ∈ VC tel que
(4.1.25)
a(u, v − u) + j(u, v) − j(u, u) ≥ l(v − u) ∀v ∈ VC

Remarque : on peut montrer que le problème (4.1.18) est équivalent au problème (4.1.25).
Cocu [Cocu 1984] a démontré que le problème de Signorini avec frottement de Coulomb
(4.1.18) admet au moins une solution si le coefficient de frottement µ est “petit”. L’unicité de
la solution est encore un problème ouvert.

4.2 Quelques méthodes pour résoudre le problème de contact

avec ou sans frottement
De très nombreuses méthodes ont été développées pour résoudre le problème de contact.
Nous présentons brièvement quatre méthodes classiques de formulation du problème de
contact i.e. la méthode de pénalisation, la méthode des multiplicateurs de Lagrange, la mé-
thode du Lagrangien augmenté et la méthode de point fixe. Ensuite, quelques algorithmes
adaptés aux différentes méthodes sont également présentés. Pour plus de détails, le lecteur
pourra se reporter à [Raous et al. 1988, Raous 1999, Wriggers 2006].

4.2.1 Méthode de pénalisation

La méthode de pénalisation consiste à approcher le problème en introduisant un
paramètre suffisamment grand permettant de définir un problème de minimisation
[Glowinski et al. 1976, Kikuchi & Song 1981, Chabrand et al. 1998].
La méthode de pénalisation remplace les lois de contact et de frottement par des lois plus
régulières qui peuvent être présentées dans les figures 4.6 et 4.7.
La solution de ce problème tend vers la solution exacte quand le paramètre de pénalisation
est suffisamment grand. Cette méthode est très simple à mettre en oeuvre. Cependant, cette
92 Chapitre 4. Application de la stratégie au problème de contact

glissement

décollement adhérence
jeu

glissement
contact
F I G U R E 4.6 – Représentation la loi de F I G U R E 4.7 – Représentation la loi de
contact de façon plus régulière Coulomb de façon plus régulière

méthode ne vérifie pas exactement les conditions de contact. De plus, le paramètre de pénali-
sation est difficile à définir : il doit être suffisamment grand pour que le problème approché
soit proche du problème réel mais s’il est trop grand, le système à résoudre devient très vite
mal conditionné.
Cette méthode est largement utilisée par les codes industriels, cf. Abaqus [Abaqus 2016]
Code_Aster [Code_Aster 2016].

4.2.2 Méthode des multiplicateurs de Lagrange

La méthode des multiplicateurs de Lagrange introduit les multiplicateurs de Lagrange
pour représenter les forces de contacts [Kikuchi & Oden 1988, Chabrand et al. 1998]. Le
problème de contact unilatéral sans frottement devient un problème de minimisation sans
contraintes et le problème de contact unilatéral avec frottement devient un problème de
pseudo-minimisation. Le problème de pseudo-minimisation peut être associé à la méthode du
Lagrangien augmenté [Alart & Curnier 1991].
Les conditions de contact sont obtenues dans ce cas de façon exacte, contrairement à la
méthode de pénalisation. Cependant, cette méthode augmente le nombre de degrés de liberté
du problème par l’introduction de nouvelles inconnues.
L’algorithme d’Uzawa peut être utilisé pour résoudre le problème de minimisation mais il
est très coûteux.

4.2.3 Méthode du Lagrangien augmenté

La méthode du Lagrangien augmenté peut être vue comme une combinai-
son de la méthode des multiplicateurs de Lagrange et de la méthode de pé-
nalisation [Alart & Curnier 1991, Simo & Laursen 1992, Heegaard & Curnier 1993,
Wriggers & Zavarise 1993, Chabrand et al. 1998]. La vitesse de convergence voire la conver-
gence de cette méthode dépend du paramètre d’augmentation. Ce type de méthode est assez
facile à intégrer dans les codes de calcul industriels, cf. Code_Aster [Code_Aster 2016].
4.2. Méthodes pour résoudre le problème de contact 93

4.2.4 Méthode de point fixe

Le problème (4.1.25) est dit implicite car le terme j(u, v) contient FCN qui dépend de la
solution u. En vue de résoudre ce problème implicite, la méthode du point fixe sur la limite de
glissement µFCN a été introduite pour transformer une inéquation quasi-variationnelle (4.1.25)
en une inéquation variationnelle plus classique [Panagiotopoulos 1985, Raous et al. 1988].
Avec la fonction de la limite de glissement g définie par le problème (4.2.1), nous résoudrons
un problème de frottement de Tresca (4.2.2) à chaque itération de point fixe.

Trouver g point fixe de l’application g 7→ µ|FCN |

(
(4.2.1)
avec u la solution du problème 4.2.2 dépendant de g
(
Trouver u ∈ VC tel que
(4.2.2)
a(u, v − u) + j(v) − j(u) ≥ l(v − u) ∀v ∈ VC
avec Z
j(v) = g|vT |dΓ
ΓC

Duvaut et Lions [Duvaut & Lions 1976] ont montré que le problème de contact statique
avec la loi de frottement de Tresca (4.2.2) peut être aussi considéré comme un problème de
minimisation : (
Trouver u ∈ VC tel que
(4.2.3)
J(u) ≤ J(v) ∀v ∈ VC
avec
1
J(v) = a(v, v) − l(v) + j(v)
2
Le problème (4.2.2) a une solution unique [Duvaut & Lions 1976]. Ainsi, par la méthode
de point fixe, le problème de Signorini avec la loi de frottement de Coulomb est résolu via des
résolutions intermédiaires d’une suite de problèmes avec la loi de frottement de Tresca.
Nous présentons quelques algorithmes adaptés aux différentes méthodes de formulation
pour résoudre le problème de contact.

4.2.5 Algorithmes de résolutions

[Link] Méthodes de projection

Les méthodes de projection consistent à projeter les conditions unilatérales sur chaque
noeud dans la zone de contact. À chaque itération de résolution, les méthodes vérifient
les conditions de contact : si le jeu devient négatif, une condition de contact est imposée.
L’itération continue jusqu’à la convergence.
La plus classique de ces méthodes est la méthode de Gauss-Seidel avec ou sans paramètre
de relaxation [Raous et al. 1988]. La méthode de point fixe peut être couplée à ce type de
méthodes.
94 Chapitre 4. Application de la stratégie au problème de contact

[Link] Méthodes des statuts

Les méthodes des statuts (en anglais “active set method”) consistent à donner a priori
les statuts de chaque noeud dans la zone de contact. Les méthodes des statuts sont souvent
utilisées en étant couplées avec la méthode des multiplicateurs de Lagrange. Les conditions
de contact seront mises à jour par un algorithme itératif. Après chaque résolution, les forces
de contact et le jeu sont analysés : si les forces de contact sont des tractions, les conditions de
contact sont supprimées et si le jeu est inférieur à 0, des conditions de contact sont ajoutées
par la suite (voir l’algorithme 4). Une comparaison de deux méthodes des statuts avec la
méthode du Lagrangien augmenté a été étudiée par [Abide et al. 2016].
Ces méthodes sont simples à mettre en oeuvre. La convergence des méthodes des statuts
est démontrée si une seule condition de contact est mise à jour par itération. Une d’entre elles
est implémentée dans Cast3M [CAST3M 2016] et sera utilisée par la suite car ces méthodes
semblent être bien adaptées à un couplage avec la méthode multi-grilles locale LDC, comme
on le verra par la suite.

[Link] Méthodes de programmation mathématique

Ce type de méthode se base sur la forme complémentaire du problème de contact avec

frottement [Lemke 1980, Klarbring 1986, Klarbring & Björkman 1988, Lebon & Raous 1992,
Chabrand et al. 1998]. Les méthodes de programmation mathématique sont efficaces et pré-
cises. Elles ont été utilisées pour résoudre des problèmes de contact en 3D. Cependant, ce
type de méthode ne converge pas toujours si la taille du système est très grande.

4.3 Stratégie de raffinement automatique du maillage adaptée au

problème de contact
Dans cette section, nous appliquons la stratégie de raffinement présentée dans le chapitre 3
au problème de contact.
L’étude de Wriggers et Scherf [Wriggers & Scherf 1998] a montré que l’estimateur ZZ fonc-
tionne correctement pour le problème de contact. Comme proposé par Nambiar et Lawrence
[Nambiar & Lawrence 1992] pour des problèmes multi-matériaux, nous estimons séparément
les erreurs sur la gaine et sur la pastille.
Comme présenté dans le section 3.2, nous générons les sous-niveaux à la première étape
de prolongation. A partir de la solution du calcul, nous estimons les erreurs séparément dans
la gaine et la pastille en utilisant l’estimateur ZZ. On choisit les éléments possédant une
erreur supérieure à l’erreur fixée par l’utilisateur pour générer un sous-niveau. Dans Cast3M,
les éléments de contact sont “noeud-facette” et nous utilisons un traitement symétrique du
contact. Nous choisissons d’élargir de façon automatisée les zones à raffiner pour que les deux
solides aient les même zones de contact et pour appliquer la méthode LDC à des maillages
réglés.
4.3. Stratégie adaptée au problème de contact 95

Algorithme 4 : Algorithme de la méthode des statuts

Initialisation : activer contact sur les noeuds dans ΓC
/*Itération*/
répéter
Résoudre le problème avec les conditions d’égalité (en déplacement) ;
Vérifier les statuts des noeuds :
si FCN > 0 alors
Supprimer la condition de contact sur ce noeud ;
fin
si uN > d alors
Ajouter la condition de contact sur ce noeud ;
fin
Mettre à jour les conditions d’égalité (en déplacement) ;
jusqu’à convergence

En utilisant la méthode des statuts, nous pouvons facilement adapter la stratégie de raf-
finement automatique de maillage proposée au problème de contact. Nous décrivons dans
l’algorithme 5 en page 96 l’utilisation de la stratégie de raffinement automatique de maillage
pour traiter un problème de contact.
Nous proposons ici de continuer à utiliser l’algorithme LDC standard, i.e. en ne re-
streignant que les déplacements, contrairement aux méthodes multi-grilles standards (voir
[Lebon 1995]). En effet, on choisit ici de résoudre complétement le problème de contact sur
toutes les grilles.
96 Chapitre 4. Application de la stratégie au problème de contact

Algorithme 5 : Algorithme LDC appliqué au problème de contact couplé avec l’estima-

teur d’erreur
Initialisation [U00 ] : résoudre le problème initial par l’algorithme 4;
Itération : calcul de [Ulk ] ;
pour k=1 à k ∗ faire
[U0k ] = [U0k−1 ];
/*Étape de prolongement*/;
pour l = 1 à l∗ faire
si k=1 alors
/* Génération des sous-niveaux et définition des conditions aux limites sur Γl ∩ Γ*/ ;
1
Application de l’estimateur d’erreur sur la solution [Ul−1 ] afin de détecter la zone à
raffiner ;
Raffinement de la zone détectée : obtention de la grille Gl ;
Définition des conditions aux limites de Dirichlet [Ul,D ] et de Neumann [Gl ] sur
Γl ∩ Γ à partir de celles définies sur Γ ;
Définition du seconde membre [Fl0 ] par discrétisation f sur Gl ;
fin
Détermination des conditions aux limites sur Γl \Γ :
k l k
[Ul,D ] = [Pl−1 ][Ul−1 ] sur Γl \Γ

Détermination des zones de contact ;

Calcul de [Ulk ] par l’algorithme 4;
fin
/*Étape de restriction*/;
pour l = l∗ − 1 à 0 faire
k
Restriction de la solution fine [Ul+1 ] sur Al :

[Ũlk ] = [Rl+1
l k
][Ul+1 ] sur Al

Calcul du résidu local sur Ål :

[Dlk ] = [Kl ][Ũlk ] − [Fl0 ] sur Ål

Correction du second membre :

[Flk ] = [Fl0 ] + χÅl [Dlk ]

Définition des zones de contact ;

Calcul de [Ulk ] par l’algorithme 4;
fin
fin
4.4. Résultats numériques 97

4.4 Résultats numériques

4.4.1 Présentation du cas test
Nous modélisons un contact entre la pastille et la gaine dans une configuration 2D axi-
symétrique. Les matériaux étudiés sont élastiques : le matériau de la pastille a un module
d’Young de 190 GPa et un coefficient de Poisson ν = 0, 3, le matériau de la gaine a un module
d’Young de 78 GPa et un coefficient de Poisson ν = 0, 34. Afin d’obtenir une inhomogénéité
du contact, une pression discontinue est appliquée sur l’extérieur de la gaine (voir figure 4.8).
Nous avons P0 = P1 = 0 et P2 = 200 MPa. L’ensemble des conditions aux limites est présenté
dans la figure 4.8.

6.4 mm

0.6 mm

4.1mm 0.57 mm

F I G U R E 4.8 – Schéma du problème

4.4.2 Illustration de la détection des zones à raffiner

Dans la section 4.3, nous avons modifié la stratégie proposée de raffinement automatique
du maillage pour l’adapter au problème de contact. Nous présentons des illustrations de la
détection des zones à raffinement avec notre cas test.
A partir de la solution initiale avec un maillage de départ de pas de 328 µm présenté dans
la figure 4.9, nous estimons séparément dans les deux corps les erreurs relatives en norme
énergie. Les erreurs estimées sont présentées dans la figure 4.10. Pour un seuil d’erreur β fixé
à 5%, les zones détectées à raffiner et les zones choisies à raffiner sont présentées dans la
figure 4.11.
Le rapport de raffinement, les éléments utilisées et les opérateurs de prolongement et de
restriction sont identiques que ceux dans la section 3.4.1.
98 Chapitre 4. Application de la stratégie au problème de contact

VONMISES
> 9.75E+05 SCAL SCAL
< 2.36E+08 > 9.77E−04 > 4.92E−03
< 2.96E−01 < 5.91E−01
2.34E+08

2.23E+08

2.12E+08 5.00E−02 5.00E−02

2.01E+08
4.50E−02 4.50E−02
1.89E+08

1.78E+08 4.00E−02 4.00E−02

1.67E+08 3.50E−02 3.50E−02

1.56E+08
3.00E−02 3.00E−02
1.45E+08

1.33E+08 2.50E−02 2.50E−02

1.22E+08
2.00E−02 2.00E−02
1.11E+08

9.99E+07 1.50E−02 1.50E−02

8.87E+07
1.00E−02 1.00E−02
7.75E+07

6.63E+07 8.20E−03 8.20E−03

5.51E+07
6.40E−03 6.40E−03
4.39E+07
4.60E−03 4.60E−03
3.27E+07

2.15E+07 2.80E−03 2.80E−03

1.03E+07
1.00E−03 1.00E−03

F I G U R E 4.9 – Exemple des

F I G U R E 4.10 – Estimation séparée des er-
contraintes sur le maillage gros-
reurs dans la gaine et la pastille
sier

Nous utilisons une solution obtenue avec maillage très fin comme solution de référence
(taille de maille environ 3,1 µm, nombre de noeuds 1653223).

4.4.3 Résultats numérique pour un cas test sans jeu ni frottement

Nous commençons par un cas test avec un jeu initial nul entre la pastille et la gaine. Nous
sommes donc dans une situation de contact affleurant.
La figure 4.12 (respectivement la figure 4.13) présente la superposition de tous les niveaux
de maillage pour un pas de maillage de départ de 328 µm (respectivement 164 µm) et
un seuil β fixé à 5%. Les figures nous montrent que la stratégie de raffinement proposée
génère bien des sous-niveaux dans les mêmes zones pour différents pas de maillage de départ.
Pour un pas de maillage local identique, la zone raffinée obtenue est identique. Il semble
donc intéressant d’utiliser cette stratégie de raffinement avec un maillage de départ assez
grossier. En visualisant la figure 4.15 qui présente les forces de contact et la figure 4.14 qui
présente la déformation pour un maillage de départ de 328 µm, nous pouvons conclure que
les sous-niveaux sont ajoutés dans la zone intermédiaire entre les zones d’adhésion et de
décollement. Dans cette zone, les erreurs sont en général importantes. Ces premiers résultats
nous confirment le bon fonctionnement de l’estimateur ZZ pour des problèmes de contact.
Nous appliquons ensuite la stratégie de raffinement à différents pas de maillages de départ
(328 µm, 164 µm, 82 µm, 41 µm) et différents seuils d’erreur β (5%, 2%, 1%). Les résultats
sont présentés dans le tableau 4.1 en cinq lignes : la première ligne indique le nombre de
sous-niveaux générés automatiquement, la deuxième ligne présente l’erreur relative entre la
solution composite à convergence et la solution de référence, la troisième ligne contient le
nombre de nœuds de chaque niveau, la quatrième ligne présente le nombre total de nœuds
et la cinquième ligne montre le temps total de calcul.
4.4. Résultats numériques 99

F I G U R E 4.11 – Choix les zones à raffiner

Seuil d’erreur
Taille de maille initiale 5% 2% 1%
nombre de sous-niveaux : 3 4 5
erreur relative : 4,08% 1,23 % 0,99%
hi nombre de noeuds par niveau : 1785/2484/5060/8090 1785/5644/12788/23030/52854 1785/6188/21038/49306/87300/202876
nombre total de noeuds : 18237 96101 368493
temps de calcul : 54s 718s 4982s
nombre de sous-niveaux : 2 3 4
erreur relative : 2,60 % 1,06% 0,78%
hi/2 nombre de noeuds par niveau : 6800/4972/8908 6800/13066/22172/52026 6800/21306/50232/86580/200732/
nombre total de noeuds : 20680 94064 365650
temps de calcul : 56s 412s 3688s
nombre de sous-niveaux : 1 2 3
erreur relative : 2,15% 1,35% 0,85%
hi/4 nombre de noeuds par niveau : 26532/8772 26532/22654/51750 26532/50778/86220/199660
nombre total de noeuds : 35304 100936 363190
temps de calcul : 106s 258s 1943s
nombre de sous-niveaux : 0 1 2
erreur relative : 2,23% 1,30% 0,85%
hi/8 nombre de noeuds par niveau : 104410 104410/51750 104410/86220/199124
nombre total de noeuds : 104410 156160 613873
temps de calcul : 300s 781s 4057s

Tableau 4.1 – Cas test de contact 2D(r,z) pour différents pas de maillages de départ et diffé-
rents seuils fixés · hi=328 µm · Sans jeu initial · Critère d’arrêt : mesure inférieure à 0,5 % de
la mesure du domaine initial.

Les résultats du tableau 4.1 montrent que la stratégie de raffinement automatique couplant
la méthode LDC et l’estimateur ZZ fonctionne très bien dans ce cas et que l’algorithme LDC
est très robuste car il marche directement sur des problèmes de contact. L’erreur composite
diminue rapidement par ajout d’un sous-niveau. De plus, on obtient le même niveau d’erreur
pour un pas de sous-niveau le plus fin équivalent. La stratégie de raffinement permet de
plus de raffiner toujours les mêmes zones quelque soit le maillage de départ (même nombre
d’éléments dans les sous-niveaux de pas de maillage équivalent). La génération des sous-
niveaux s’arrête automatiquement. Le critère d’arrêt géométrique avec 0,5% de la mesure du
100 Chapitre 4. Application de la stratégie au problème de contact

F I G U R E 4.12 – Superposition des ni- F I G U R E 4.13 – Superposition des ni-

veaux de maillage pour un maillage initial veaux de maillage pour un maillage initial
de pas 328 µm et un seuil fixé à 5% · Sans de pas 164 µm et un seuil fixé à 5% · Sans
jeu initial. jeu initial.

domaine initial semble être assez générique. A la fin des calculs, les erreurs entre les solutions
composites et la solution de référence respectent les seuils demandés.
Les résultats montrent l’efficacité de la méthode LDC. Pour un seuil fixé à 5%, l’ajout de
sous-niveaux permet de diminuer le temps de calcul tout en ayant le même niveau de précision
qu’en utilisant d’un maillage très fin uniforme. Pour les seuils fixés à 2% et 1%, il semble qu’il
y ait à faire un compromis entre le nombre de sous-niveaux et temps de calcul. Les résultats
les plus optimaux en temps de calcul sont obtenus pour un maillage de départ avec un pas
de hi/4. Afin de comprendre cette tendance, nous avons superposé les sous-niveaux générés
pour le maillage de départ de hi et le seuil fixé à 1% sur la figure 4.16. Nous voyons que
pour atteindre ce niveau de précision, nous avons besoin de raffiner presque tout le domaine.
Le maillage de départ de hi est donc trop grossier pour atteindre rapidement le niveau de
précision demandé. A partir du troisième sous-niveau dont la taille de maille correspond au
maillage de départ de hi/4, les sous-niveaux commencent à être ajoutés localement, et c’est
là que la méthode LDC devient efficace.
4.4. Résultats numériques 101

F I G U R E 4.14 – Déformation pour un F I G U R E 4.15 – Forces de contact pour un

maillage initial de pas 328 µm (facteur maillage initial de pas 328 µm · Sans jeu
d’amplification = 50) · Sans jeu initial. initial.

Optimisation des temps de calcul

Nous proposons de regarder la répartition des temps de calcul en trois parties pour les
calculs avec les seuils β = 1% et 2% :
• Temps de résolution initiale du maillage grossier
• Temps de génération des sous-niveaux qui correspond au temps de la première étape
de prolongement
• Temps des ∧-cycles depuis la première restriction jusqu’à la convergence du calcul
Les répartitions des temps de calcul sont présentées sur les figures 4.17 et 4.18.
Les figures 4.17 et 4.18 montrent que les ∧-cycles pour les problèmes de contact semblent
être assez coûteux. En effet, la non-linéarité de contact doit être résolue itérativement à
chaque appel du solveur. Comme le nombre de ∧-cycles est proportionnel au nombre de
sous-niveaux, le temps de calcul des ∧-cycles sera proportionnel au nombre de sous-niveaux,
on a donc intérêt à limiter le nombre de sous-niveaux. A contrario, le temps de calcul sur le
maillage initial augmente avec la diminution du pas du maillage. Il y a donc un optimum
« nombre de sous niveaux/pas du maillage initial » à trouver. Ici l’optimum semble obtenu
pour un maillage initial de pas hi dans le cas d’une erreur voulue de 5% et de pas hi/4 pour
une erreur voulue de 2% ou 1%.
En vue de diminuer le temps de calcul des ∧-cycles, nous proposons d’utiliser une tech-
nique permettant d’initialiser la zone de contact actif : au lieu d’initialiser les contacts actifs
sur toute la zone possible de contact, on active les statuts de contact uniquement sur la zone de
102 Chapitre 4. Application de la stratégie au problème de contact

F I G U R E 4.16 – Superposition des niveaux de maillage pour un maillage initial de pas 328 µm
et un seuil fixé à 1% · Sans jeu initial.

800

700 calcul -cycle

v
calcul première génération
600 calcul grossier
44%

500
Temps (s)

75%
400
19%
300 57%

200 35%

37%
100 25% 36% 48%

0.095% 6.8% 17%

0
hi hi/2 hi/4 hi/8
Taille du maillage grossier

F I G U R E 4.17 – Répartition des temps de calcul pour β = 2%

contact obtenue après résolution du problème immédiatement plus grossier dans la première
étape de prolongement. Pendant les ∧-cycles, l’initialisation des statuts de contact sera définie
à partir de la solution précédente de même niveau pour éviter une extrapolation de la zone de
contact et optimiser la convergence de la méthode des statuts à convergence de l’algorithme
LDC. Dans la suite, cette technique sera appelée “maillage de contact actif” (MCA). Nous
4.4. Résultats numériques 103

calcul -cycle
v
5000 calcul première génération
calcul grossier

4000
Temps (s)

73%
3000
61%
60%

2000
27%

1000 32%
27% 39% 71%

0.014% 0.81% 2.3% 7%

0
hi hi/2 hi/4 hi/8
Taille du maillage grossier

F I G U R E 4.18 – Répartition des temps de calcul pour β = 1%

testons cette technique sur différents pas de maillages de départ (hi=328 µm, hi/2=164 µm,
hi/4=82 µm, hi/8=41 µm) et différents seuils d’erreur β (5%, 2%, 1%). Les résultats sont
présentés dans le tableau 4.2.
Seuil d’erreur
Taille de maille initiale 5% 2% 1%
nombre de sous-niveaux : 3 4 5
erreur relative : 4,08% 1,23 % 0,99%
hi nombre de noeuds par niveau : 1785/2484/5060/8090 1785/5644/12788/23030/52854 1785/6188/21038/49306/87300/202876
nombre total de noeuds : 18237 96101 368493
temps de calcul : 44s 417s 2764s
nombre de sous-niveaux : 2 3 4
erreur relative : 2,60 % 1,06% 0,78%
hi/2 nombre de noeuds par niveau : 6800/4972/8908 6800/13066/22172/52026 6800/21306/50232/86580/200732/
nombre total de noeuds : 20680 94064 365650
temps de calcul : 51s 320s 1627s
nombre de sous-niveaux : 1 2 3
erreur relative : 2,15% 1,35% 0,85%
hi/4 nombre de noeuds par niveau : 26532/8772 26532/22654/51750 26532/50778/86220/199660
nombre total de noeuds : 35304 100936 363190
temps de calcul : 113s 207s 986s
nombre de sous-niveaux : 0 1 2
erreur relative : 2,23% 1,30% 0,85%
hi/8 nombre de noeuds par niveau : 104410 104410/51750 104410/86220/199124
noeud total : 104410 156160 389754
temps de calcul : 300s 780s 1694s

Tableau 4.2 – Cas test de contact 2D(r,z) sans frottement pour différents pas de maillages de
départ et différents seuils fixés · hi=328 µm · Sans jeu initial · Initialisation des statuts de
contact par la technique MCA · Critère d’arrêt : mesure inférieure à 0,5% de la mesure du
domaine initial.

Les résultats présentés dans le tableau 4.2 sont similaires aux résultats présentés dans le
tableau 4.1 quant au nombre de sous-niveaux, à l’erreur et au nombre de noeuds. La solu-
104 Chapitre 4. Application de la stratégie au problème de contact

800

calcul -cycle
v
700
calcul première génération
calcul grossier
600

54%
500
Temps (s)

400

300 8.2%
72%
59%
200
32%
38%
100 32% 46%
28%
0.16% 8.7% 22%
0
hi hi/2 hi/4 hi/8
Taille du maillage grossier

F I G U R E 4.19 – Répartition des temps de calcul pour β = 2% pour la technique MCA

4000
calcul -cycle
v
3500 calcul première génération
calcul grossier
3000

2500
Temps (s)

79%
2000

1500
54% 48%
1000
34%
500 35%
21% 44% 62%
0.019% 1.8% 4.5% 16%
0
hi hi/2 hi/4 hi/8
Taille du maillage grossier

F I G U R E 4.20 – Répartition des temps de calcul pour β = 1% pour la technique MCA

tion du problème est donc bien inchangée. Cette technique par ailleurs a beaucoup diminué
(jusqu’à un facteur 2) le temps de calcul pour les seuils fixés à 2% et 1%. Comme attendu, il
est intéressant d’utiliser la technique MCA quand il y a plusieurs sous-niveaux afin de trans-
mettre des informations de contact entre les ∧-cycles. Nous allons conserver cette technique
d’initialisation des statuts de contact dans la suite.
Les répartitions des temps de calcul avec la technique MCA sont présentées sur les fi-
gures 4.19 et 4.20. La tendance n’a pas trop été modifiée : il y a toujours un minimum à
4.4. Résultats numériques 105

trouver entre le maillage de départ et la localisation des sous-niveaux.

4.4.4 Résultats numériques pour un cas test avec jeu sans frottement
Lors de l’irradiation du réacteur, il y a un jeu initial entre la pastille et la gaine d’environ
85 µm. Ce jeu se ferme à la fin du premier cycle de vie en réaction (voir section 1.4). Nous
souhaitons tester la stratégie de raffinement à un cas test avec un jeu initial. Ce jeu est choisi
petit afin de représenter l’instant de contact. Ici on prend un jeu initial de 2 µm.
La figure 4.21 présente la superposition des niveaux de maillage pour un maillage initial
de pas hi=328 µm et un seuil fixé β à 5% pour le cas test avec jeu initial égal à 2 µm. Les
sous-niveaux se localisent bien autour de la zone de changement de statut de contact. La
stratégie proposée peut donc générer des sous-niveaux dans la zone où l’erreur est importante
pour un cas test avec jeu. La figure 4.22 montre les forces de contact pour ce cas test. Les
forces présentent bien seulement des composantes normales.

F I G U R E 4.21 – Superposition des ni- F I G U R E 4.22 – Forces de contact pour un

veaux de maillage pour un maillage initial maillage initial de pas 328 µm pour le cas
de pas 328 µm et un seuil fixé à 5% · Jeu test · Jeu initial égal à 2 µm.
initial égal à 2 µm.

Des résultats plus détaillés pour différents pas de maillages de départ (328 µm, 164 µm,
82 µm, 41 µm) et différents seuils d’erreur β (5%, 2%, 1%) sont présentés dans le tableau 4.3.
Les résultats présentés dans le tableau 4.3 montrent que la stratégie de raffinement fonc-
tionne très bien pour un cas test de contact avec jeu initial. La stratégie permet d’ajouter des
sous-niveaux dans la zone où les erreurs sont les plus importantes et d’arrêter de générer les
sous-niveaux quand les solutions sont suffisamment précises en regard des seuils fixés. À la
fin des calculs, les erreurs entre les solutions composites et la solution de référence respectent
106 Chapitre 4. Application de la stratégie au problème de contact

Seuil d’erreur
Taille de maille initiale 5% 2% 1%
nombre de sous-niveaux : 3 4 5
erreur relative : 3,14% 1,77% 1,07%
hi nombre de noeuds par niveau : 1785/2484/5500/13662 1785/5508/13152/22842/55614 1785/6800/21306/49500/86580/199124
nombre total de noeuds : 23431 98901 365095
temps de calcul : 21s 315s 1754s
nombre de sous-niveaux : 2 3 4
erreur relative : 2,03 % 1,01% 0,6%
hi/2 nombre de noeuds par niveau : 6800/5500/14076 6800/13700/22842/55338 6800/21038/50948/85500/198990
nombre total de noeuds : 26376 98680 363276
temps de calcul : 32s 209s 1669s
nombre de sous-niveaux : 1 2 3
erreur relative : 1,97% 1,36% 0,83%
hi/4 nombre de noeuds par niveau : 26532/13662 26532/22654/55062 26532/52032/85500/196444
nombre total de noeuds : 40194 104248 360508
temps de calcul : 82s 190s 912s
nombre de sous-niveaux : 0 1 2
erreur relative : 2,12% 1,47% 0,45%
hi/8 nombre de noeuds par niveau : 104410 104410/54264 104410/85140/196444
nombre total de noeuds : 104410 158674 385994
temps de calcul : 300s 754s 2350s

Tableau 4.3 – Cas test contact 2D(r,z) sans frottement pour différents pas de maillages de
départ et différents seuils fixés · hi=328 µm · Jeu initial = 2µm · Initialisation des statuts de
contact par la technique MCA · Critère d’arrêt : mesure inférieure à 0,5% de la mesure du
domaine initial.

bien les seuils demandés. Ici encore, les zones de raffinement obtenues sont les mêmes pour
un pas de maillage local identique. Il existe encore un optimum entre maillage de départ et
localisation des sous-niveaux (et donc nombre de ∧-cycles). Ici aussi, le maillage de départ
hi/4 semble optimal pour les seuil fixés à 2% et 1%.

4.4.5 Résultats numériques pour un cas test avec jeu avec frottement
Nous voulons à présent tester notre stratégie de raffinement automatique de maillage pour
un problème de contact avec frottement de Coulomb. Le coefficient de frottement est égal à
0,2 (µ = 0.2).
Les sous-niveaux générés pour un maillage initial de pas 328 µm et un seuil fixé à 5% sont
superposés dans la figure 4.23. La stratégie génère toujours bien les sous-niveaux autour de
la zone de changement de statuts de contact.
La figure 4.24 représente les forces de contact avec frottement. Contrairement aux cas
tests sans frottement, les forces de contact ont des composantes normales et des composantes
tangentielles. Un glissement peut être observé dans la zone de contact.
Le tableau 4.4 présente des résultats de calcul pour différents pas de maillages de départ
(328 µm, 164 µm, 82 µm, 41 µm) et différents seuils d’erreur β (5%, 2%, 1%).
Les résultats montrent que la stratégie de raffinement nous permet d’obtenir de bonnes
performances également sur les problèmes de contact avec frottement. La stratégie génère
automatiquement des sous-niveaux dans les zones d’intérêt et s’arrête automatiquement de
générer des sous-niveaux lorsque la précision demandés est atteinte. Ainsi, nous atteignons le
même niveau d’erreur pour des tailles de maille de niveaux les plus fins égales. Au niveau du
4.4. Résultats numériques 107

temps de calcul, il y a là encore un optimum qui dépend du pas du maillage initial et du seuil
d’erreur. Par exemple, il semble que le maillage de départ avec un pas de hi/2 soit optimal
pour les seuils 5% et 2% alors qu’un pas de hi/4 est optimal pour un seuil de 1%.

F I G U R E 4.23 – Superposition des niveaux de maillage pour un maillage initial de pas 328 µm
et un seuil fixé à 5% · Jeu initial égal à 2 µm · Coefficient de frottement égal à 0,2.

F I G U R E 4.24 – Forces de contact avec frottement pour un maillage initial de pas 328 µm ·
Jeu initial égal à 2 µm · Coefficient de frottement égal à 0,2. Figure de droite zoom de la zone
en rouge
108 Chapitre 4. Application de la stratégie au problème de contact

Seuil d’erreur
Taille de maille initiale 5% 2% 1%
nombre de sous-niveaux : 3 4 5
erreur relative : 6,04% 1,27% 0.84%
hi nombre de noeuds par niveau : 1785/2484/4972/11022 1785/5508/12604/21988/51000 1785/6800/21038/48048/83952/194142
nombre total de noeuds : 20263 92885 355765
temps de calcul : 90s 765s 6477s
nombre de sous-niveaux : 2 3 4
erreur relative : 2,79 % 1,06% 0,51%
hi/2 nombre de noeuds par niveau : 6800/5500/14076 6800/12878/21804/49714 6800/21038/48958/83838/194040
nombre total de noeuds : 22706 91196 354674
temps de calcul : 81s 506s 4297s
nombre de sous-niveaux : 1 2 3
erreur relative : 2,37% 1,31% 0.89%
hi/4 nombre de noeuds par niveau : 26532/13662 26532/21804/49446 26532/49496/82544/192456
nombre total de noeuds : 37554 97782 351028
temps de calcul : 170s 596s 4152s
nombre de sous-niveaux : 0 1 2
erreur relative : 2,15% 1,41% 0,42%
hi/8 nombre de noeuds par niveau : 104410 104410/49446 104410/82544/190474
nombre total de noeuds : 104410 153856 377428
temps de calcul : 591s 1422s 5788s

Tableau 4.4 – Cas test contact 2D(r,z) avec frottement pour différents pas de maillages de
départ et différents seuils fixés · hi=328 µm · Jeu initial = 2µm · Coefficient de frottement
égal à 0,2 · Initialisation des statuts de contact par la technique MCA · Critère d’arrêt : mesure
inférieure à 0,5% de la mesure du domaine initial.

Calculs avec différents coefficients de frottement

Pour finir, nous allons tester la stratégie proposée pour différents coefficients de frottement.
Nous reprenons le cas test précédent avec des coefficients de frottement µ = 0, 6 ou 1, 0. Les
forces de contact dans ces deux cas, pour un maillage initial de pas 328 µm, sont présentées
dans les figures 4.25 et 4.26. La direction des forces montre que les forces tangentielles sont
plus petites que µFcN . Il n’y a pas de glissement entre les deux objets.
Nous présentons les calculs pour différents pas de maillages de départ (328 µm, 164 µm,
82 µm, 41 µm) et différents coefficients de frottement µ (0,2 ; 0,6 ; 1,0) dans le tableau 4.5.
Le seuil β est fixé à 5%.
Les résultats présentés dans le tableau 4.5 montrent que la stratégie proposée de raffi-
nement fonctionne bien pour différents coefficients de frottement. La stratégie semble être
très générique. Les solutions ont les mêmes nombres de sous-niveaux et les mêmes niveaux
d’erreur pour des calculs de même pas de maillage de départ. Les erreurs sont inférieures aux
seuils demandés sauf pour le pas de maillage de 328 µm. Il semble que la stratégie de raffine-
ment détecte presque les mêmes zones pour les coefficients de frottement testés. Cependant,
il semble que le temps de calcul augmente quand le coefficient de frottement augmente. Ceci
semble dû à la convergence de la méthode des statuts (voir la ligne correspondant à un pas
de maillage de départ de hi/8 pour lequel l’algorithme LDC n’a pas été appliqué). Comme
pour tous les cas tests précédents, il existe aussi un compromis entre la taille de maillage de
départ et le nombre de sous-niveaux.
4.4. Résultats numériques 109

F I G U R E 4.25 – Forces de contact avec F I G U R E 4.26 – Forces de contact avec

frottement pour un maillage initial de pas frottement pour un maillage initial de pas
328 µm · Jeu initial égal à 2µm · Coeffi- 328 µm · Jeu initial égal à 2µm · Coeffi-
cient de frottement égal à 0,6. cient de frottement égal à 1,0.

Coefficient de frottement
Taille de maille initiale 0,2 0,6 1,0
nombre de sous-niveaux : 3 3 3
erreur relative : 6,04% 6,93% 6,95%
hi nombre de noeuds par niveau : 1785/2484/4972/11022 1785/2278/4746/11072 1785/2208/4522/10266
nombre total de noeuds : 20263 19881 18781
temps de calcul : 90s 68s 108s
nombre de sous-niveaux : 2 2 2
erreur relative : 2,79 % 3,37% 3,29%
hi/2 nombre de noeuds par niveau : 6800/5500/14076 6800/4662/11072 6800/4522/10266
nombre total de noeuds : 22706 22534 21588
temps de calcul : 81s 90s 97s
nombre de sous-niveaux : 1 1 1
erreur relative : 2,37% 2,98% 3,08%
hi/4 nombre de noeuds par niveau : 26532/13662 26532/10944 26532/10266
nombre total de noeuds : 37554 37476 36798
temps de calcul : 170s 210s 234s
nombre de sous-niveaux : 0 0 0
erreur relative : 2,15% 2,22% 2,27%
hi/8 nombre de noeuds par niveau : 104410 104410 104410
nombre total de noeuds : 104410 104410 104410
temps de calcul : 591s 800s 981s

Tableau 4.5 – Cas test contact 2D(r,z) avec frottement pour différents pas de maillages de
départ et différents coefficients de frottement · Seuil d’erreur β = 5% · hi=328 µm · Jeu initial
= 2µm · Initialisation des statuts de contact par la technique MCA · Critère d’arrêt : mesure
inférieure à 0,5% de la mesure du domaine initial.
110 Chapitre 4. Application de la stratégie au problème de contact

4.4.6 Contact avec frontière courbe

Dans cette section, nous nous intéressons à un cas de contact avec frontière courbe. Nous
choisissons un cas test en 2D vérifiant les hypothèses des déformations planes. Ce cas test
utilise les géométries de la pastille et de la gaine dans le plan (r, θ), cf. figure 4.28. Il s’agit
donc d’une portion du disque en face d’une portion d’anneau.
Le contact pastille-gaine est provoqué par une dilatation thermique. Afin d’éviter le calcul
thermique, nous imposons un champ de température dans les 2 corps. Les déformations
associées à ce champ de température sont calculées par ǫ = αT ∆T avec αT le coefficient
de dilatation thermique. Le champ de forces nodales résultant de l’intégration du champ de
contraintes est alors appliqué comme terme source dans le calcul mécanique.
Le coefficient de dilatation thermique de la pastille est αT = 1.0 × 10−5 et αT = 6.0 × 10−6
pour la gaine. Le champ de température est radial dans la pastille et constant dans la gaine,
voir figure 4.27. Il s’écrit sous la forme :

∆T = 7.0 × 107 × (R2 − (x2 + y 2 )) + 561. − 293. dans la pastille

∆T = 561. − 293. dans la gaine

Le jeu initial entre la pastille et la gaine est 2 µm. La figure 4.28 présente l’ensemble des
conditions aux limites.
VAL − ISO
> 2.68E+02
< 1.44E+03

1.44E+03

1.38E+03

1.32E+03

1.27E+03

1.21E+03

1.16E+03
ue
1.10E+03
étriq
1.04E+03 Sym
9.87E+02

9.31E+02

8.75E+02
ue
étriq
Sym
8.19E+02

7.63E+02

7.07E+02

6.51E+02

5.95E+02

5.39E+02

4.83E+02

4.27E+02

3.71E+02

3.15E+02 4.1mm 0.57 mm

F I G U R E 4.27 – Champ de température imposé F I G U R E 4.28 – Schéma du problème 2D

2D déformation plane déformation plane

[Link] Raffinement LDC avec maillage hiérarchique

Dans cette section, nous choisissons d’utiliser l’algorithme LDC avec raffinement hiérar-
chique, voir section 3.1. Le raffinement hiérarchique utilisé consiste à ajouter des noeuds au
milieu des segments du niveau immédiatement plus grossier et à construire les mailles fines
en s’appuyant sur ces noeuds, voir figure 4.29. On voit ainsi que les noeuds frontières ajoutés
sur le sous-niveau sont situés sur la corde du cercle et non directement sur l’arc de cercle.
Après application de l’algorithme LDC, nous vérifions les forces de contact sur chaque
niveau. La force de contact du maillage grossier est présentée dans la figure 4.30.
4.4. Résultats numériques 111

noeud de niveau fin

noeud de niveau grossier

F I G U R E 4.29 – Raffinement des mailles rectangulaires de façon hiérarchique

F I G U R E 4.30 – Forces de contact sur le maillage grossier.

On constate cependant que la force de contact du deuxième sous-niveau est très différente
de celle du maillage grossier, cf. figure 4.31 et présente des discontinuités de contact, voire
du décollement sur les noeuds fins rajoutés sur la corde.
A titre de comparaison, nous proposons de calculer le champ de force de contact obtenu

F I G U R E 4.31 – Forces de contact sur le deuxième sous-niveau. Raffinement hiérarchique.

sur un maillage initial construit de façon hiérarchique à partir du maillage grossier 4.30 et
donc équivalent localement au maillage de deuxième sous-niveau obtenu avec LDC et sans
utiliser de procédure de raffinement automatique. Les forces de contact sont présentées sur
la figure 4.32.
En comparant la figure 4.32 à la figure 4.31, on peut constater le bon fonctionnement de
la méthode LDC. En effet, les champs de force obtenus localement par la méthode LDC sont
équivalents à ceux obtenus sans LDC sur un maillage global construit de manière équivalente.
En analysant le champ de force sur le maillage déformé, cf. figure 4.33, on constate que les
112 Chapitre 4. Application de la stratégie au problème de contact

F I G U R E 4.32 – Forces de contact sur un maillage construit de façon hiérarchique à partir du

maillage grossier 4.30.

conditions aux limites sur l’axe de fissuration de la pastille (y = 0) induisent un déplacement

de celle-ci selon l’axe θ. Le contact se produit alors entre un noeud et le segment opposé sur la
gaine. Le raffinement hiérarchique avec noeuds sur la corde induit des forces normales selon
la même direction pour un segment donné. Ainsi les nouveaux noeuds générés sur la corde
du cercle ne peuvent plus rentrer en contact avec la gaine (les deux frontières de contact ne
s’imbriquent plus).

F I G U R E 4.33 – Forces de contact sur un maillage déformé 4.32

La figure 4.34 représente de manière illustrative ce phénomène. Le noeud "a" sur l’arc de
cercle de la pastille va entrer en contact avec le segment "1" et le noeud “b” avec le segment
“2”. La géométrie du maillage des 2 corps peut alors empêcher le noeud “c” de rentrer en
contact avec la gaine.

a 2
c

noeud de niveau fin

noeud de niveau grossier

F I G U R E 4.34 – Schéma illustrant le contact dans le cas d’un raffinement hiérarchique

4.4. Résultats numériques 113

[Link] Raffinement avec maillage quasi-hiérarchique

Afin de pallier le défaut du raffinement hiérarchique sur un maillage courbe de contact,

nous proposons de déplacer les noeuds ajoutés de façon hiérarchique dans la zone de contact
courbe de la corde vers l’arc de cercle, voir figure 4.35. Les maillages ainsi générés sont donc
“quasi-hiérarchiques”. La frontière de contact est donc mieux représentée à chaque niveau de
raffinement.
Nous nous limitons à utiliser ce raffinement modifié uniquement dans les zones de contact,
car avec la stratégie de détection de zone d’intérêt proposée, la solution va être recalculée
complètement dans ces zones. Sur les frontières extérieures (avec CL de Dirichlet ou Neu-
mann) ce type de raffinement poserait la question de l’extrapolation des déplacements et au
vu des résultats du chapitre 3 3, sur ces frontières le raffinement hiérarchique n’engendre pas
des erreurs du premier ordre.

noeud de niveau fin

noeud de niveau grossier

F I G U R E 4.35 – Raffinement sur l’arc de cycle de façon quai-hiérarchique

Les forces de contact obtenues en utilisant le raffinement quasi-hiérarchique sur les zones
en contact sont données sur la figure 4.36.

(a) Maillage grossier (b) Premier sous-niveau

F I G U R E 4.36 – Forces de contact après algorithme LDC et raffinement quasi-hiérarchique

sur les frontières de contact.

Les résultats montrent que le déplacement des noeuds ajoutés dans la zone de contact
vers l’arc de cercle améliore la solution. Les forces de contact semblent cohérentes sur les
sous-niveaux.
Nous avons effectué des calculs avec différents pas de maillages de départ et différents
seuils. La stratégie quasi-hiérarchique proposée permet d’obtenir des erreurs après cycles LDC
qui sont inférieures aux seuils comme le montre le tableau 4.6.
114 Chapitre 4. Application de la stratégie au problème de contact

(c) Deuxième sous-niveau (d) Troisième sous-niveau

F I G U R E 4.36 – Forces de contact après algorithme LDC et raffinement quasi-hiérarchique sur

les frontières de contact (suite)

En conclusion, le raffinement de maillage en cas de frontière de contact non-linéaire par

morceau n’est pas trivial. Le raffinement hiérarchique ne permet pas de trouver les forces de
contact souhaitées.
Une alternative quasi-hiérarchique a été proposée ici qui permet de converger vers la solution
de référence. La stratégie LDC-ZZ de raffinement automatique de maillage fonctionne encore
très bien même si elle semble conduire à trop raffiner dans ce cas. Le couplage de LDC avec
un autre estimateur plus fiable permettrait certainement de pallier ce sur-raffinement.

Bilan
Dans ce chapitre, nous avons présenté le problème de contact de Signorini avec et sans
frottement ainsi que les méthodes classiques pour le résoudre. En utilisant la méthode des
statuts, nous avons adapté la stratégie de raffinement automatique de maillage couplant la
méthode LDC et l’estimateur ZZ au problème de contact. Les résultats présentés pour des
problèmes de contact unilatéral avec ou sans jeu initial et un problème de frottement avec
jeu initial montrent que la méthode LDC et la stratégie de raffinement sont très génériques et
très efficaces. Nous disposons à présent d’un outil permettant de raffiner de façon optimale
un maillage initial d’un problème de contact en respectant un seuil d’erreur relative donnée
par l’utilisateur.
4.4. Résultats numériques 115

Seuil d’erreur
Taille de maille initiale 5% 2% 1%
nombre de sous-niveaux : 0 1 3
hi erreur relative : 1,22% 0,80% 0.47%
nombre de noeuds par niveau : 455 455/650 455/900/2254/4462
nombre de sous-niveaux : 0 0 2
hi/2 erreur relative : 0,62 % 0,62% 0,38%
nombre de noeuds par niveau : 1632 1632 1632/2162/4278
nombre de sous-niveaux : 0 0 1
hi/4 erreur relative : 0,31% 0,31% 0.26%
nombre de noeuds par niveau : 6298 6298 6298/4278
nombre de sous-niveaux : 0 0 0
hi/8 erreur relative : 0,17% 0,17% 0,17%
nombre de noeuds par niveau : 24380 24380 24380

Tableau 4.6 – Cas test contact 2D(r,θ) pour différents pas de maillages de départ et différents
seuils fixés · Jeu initial = 2µm · Raffinement avec maillage quasi-hiérarchique · Critère d’arrêt :
mesure inférieure à 0,5% de la mesure du domaine initial.
Conclusion et Perspectives

Dans l’objectif de proposer une méthode de raffinement de maillage permettant de simuler

précisément le phénomène d’Interaction mécanique Pastille-Gaine (IPG), nous avons déve-
loppé une stratégie de raffinement local automatique basée sur la combinaison de la méthode
multi-grilles locale “Local Defect Correction” (LDC) et de l’estimateur d’erreur a posteriori
de Zienkiewicz et Zhu (ZZ). L’intérêt des méthodes multi-grilles locales est de permettre
d’atteindre des tailles locales de maille très fines tout en conservant des maillages structurés
uniformes avec peu de degrés de liberté. Nous avons choisi d’utiliser la méthode LDC car
elle s’applique facilement dans le cadre d’une résolution par éléments finis et parce qu’elle
n’est pas intrusive vis-à-vis des codes de calcul existants. Afin de détecter les zones d’intérêt à
raffiner finement, l’estimateur ZZ a été choisi en raison de sa simplicité de mise en œuvre, de
sa rapidité et de son efficacité démontrée dans la littérature.
Nous avons également proposé plusieurs critères d’arrêt afin de permettre l’arrêt de la
génération de sous-niveaux lorsque le problème présente des singularités (de contraintes par
exemple). Plusieurs cas test bi-dimensionnels (axisymétriques, déformations planes) s’intéres-
sant à la réponse élastique de la gaine sont étudiés. Le contact entre la gaine et la pastille est
alors représenté par des chargements discontinus sur la paroi interne de la gaine.
Les résultats numériques basés sur l’erreur relative en norme énergie et l’erreur absolue en
norme infinie sur la contrainte équivalente de von Mises montrent que la stratégie LDC-ZZ
fonctionne très bien. Les sous-niveaux sont ajoutés automatiquement autour des zones de dis-
continuité du chargement. Les erreurs entre la solution LDC composite (rassemblant tous les
sous-niveaux) et la solution de référence sont généralement inférieures aux seuils demandés.
Parmi les différents critères d’arrêt, le critère d’arrêt basé sur un ratio de mesure de surface
entre la zone à raffiner et le domaine initial semble être générique et efficace.
Ensuite, nous avons adapté la stratégie de raffinement automatique de maillage couplant
la méthode LDC et l’estimateur ZZ à la problématique de contact avec ou sans frottement.
En utilisant la méthode des statuts pour résoudre le problème de contact, la stratégie de
raffinement proposée considérant le solveur comme une “boîte noire” s’adapte facilement
au problème de contact multi-corps. L’efficacité de la méthode LDC et de la stratégie de
raffinement proposée est montrée sur des résultats numériques de résolution d’un problème
de contact unilatéral avec ou sans jeu initial et sur un problème de frottement avec jeu initial.
La stratégie de raffinement automatique de maillage couplant la méthode LDC, l’estima-
teur ZZ et le critère d’arrêt de génération des sous-niveaux semble donc très générique et
ouvre la voie à de nouvelles applications.
118 CONCLUSION ET PERSPECTIVES

Une perspective à court terme de ce travail est de compléter l’étude sur le problème de
contact. Les résultats numériques montrent qu’il y a un compromis à trouver entre maillage
de départ et localisation des sous-niveaux pour que le temps de calcul soit optimal. Un critère
pourrait être proposé afin d’éviter de faire des sous-niveaux dans l’algorithme LDC lorsque la
zone de raffinement est trop large. Ces sous-niveaux seraient alors remplacés par un maillage
initial directement plus raffiné. Ensuite, il faudrait étendre la stratégie proposée à des pro-
blèmes comportant des matériaux non linéaires et des problèmes d’évolutions temporelles.
Les zones de raffinement évoluant alors dans le temps, il faudra s’intéresser au problème de
transmission d’informations d’un pas de temps sur l’autre entre sous-grilles et au raffinement
local spatio-temporel.
L’objectif in fine étant de disposer d’un solveur multi-grilles local pour les études concer-
nant l’interaction mécanique pastille-gaine, la perspective à plus long terme est d’intégrer
cette stratégie de raffinement automatique dans la plate-forme logicielle PLEIADES afin de
pouvoir tester l’efficacité et les performances de l’approche proposée sur des configurations
industrielles.
Bibliographie

[Abaqus 2016] Abaqus. Abaqus Analysis User’s Manual, 2016.

[Abide et al. 2016] S. Abide, M. Barboteu et D. Danan. Analysis of two active set type methods
to solve unilateral contact problems. Applied Mathematics and Computation, vol. 284,
pages 286–307, 2016.
[Ainsworth & Oden 1997] M. Ainsworth et J.T. Oden. A posteriori error estimation in finite
element analysis. Computer Methods in Applied Mechanics and Engineering, vol. 142,
pages 1–88, 1997.
[Ainsworth & Oden 2011] M. Ainsworth et J.T. Oden. A posteriori error estimation in finite
element analysis. John Wiley & Sons, 2011.
[Ainsworth et al. 1989] M. Ainsworth, J.Z. Zhu, A.W. Craig et O.C. Zienkiewicz. Analysis of
the Zienkiewicz–Zhu a-posteriori error estimator in the finite element method. Interna-
tional Journal for Numerical Methods in Engineering, vol. 28, no. 9, pages 2161–2174,
1989.
[Alart & Curnier 1991] P. Alart et A. Curnier. A mixed formulation for frictional contact pro-
blems prone to Newton like solution methods. Computer Methods in Applied Mechanics
and Engineering, vol. 92, no. 3, pages 353–375, 1991.
[Alberman et al. 1997] A. Alberman, M. Roche, P. Couffin, S. Bendotti, D.J. Moulin et J.L.
Boutfroy. Technique for power ramp tests in the ISABELLE 1 loop of the OSIRIS reactor.
Nuclear Engineering and Design, vol. 168, pages 293–303, 1997.
[Anthonissen et al. 1998] M.J.H. Anthonissen, B. van ’t Hof et A.A. Reusken. A Finite Volume
Scheme for Solving Elliptic Boundary Value Problems on Composite Grids. Computing,
vol. 61, pages 285–305, 1998.
[Anthonissen et al. 2005] M.J.H. Anthonissen, B.A.V. Bennet et M.D. Smooke. An adaptive
multilevel local defect correction technique with application to combustion. Combustion
Theory and Modelling, vol. 9, pages 273–299, 2005.
[Askes & Rodriguez-Ferran 2001] H. Askes et A. Rodriguez-Ferran. A combined rh-adaptive
scheme based on domain subdivision. Formulation and linear examples. International
Journal for Numerical Methods in Engineering, vol. 51, no. 3, pages 253–273, 2001.
[Babuška & Miller 1987] I. Babuška et A. Miller. A feedback finite element method with a
posteriori error estimation : Part I. The finite element method and some basic properties
of the a posteriori error estimator. Computer Methods in Applied Mechanics and
Engineering, vol. 61, no. 1, pages 1–40, 1987.
120 Bibliographie

[Babuška & Rheinboldt 1978a] I. Babuška et W. Rheinboldt. Error Estimates for Adaptive
Finite Element Computations. SIAM Journal on Numerical Analysis, vol. 15, no. 4,
pages 736–754, 1978.
[Babuška & Rheinboldt 1978b] I. Babuška et W.C. Rheinboldt. A-posteriori error estimates for
the finite element method. International Journal for Numerical Methods in Engineering,
vol. 12, pages 1597–1615, 1978.
[Babuška & Rheinboldt 1981] I. Babuška et W. Rheinboldt. A Posteriori Error Analysis of
Finite Element Solutions for One-Dimensional Problems. SIAM Journal on Numerical
Analysis, vol. 18, no. 3, pages 565–589, 1981.
[Babuška & Strouboulis 2001] I. Babuška et T. Strouboulis. The finite element method and
its reliability. Oxford University Press, 2001.
[Babuška & Suri 1987a] I. Babuška et M. Suri. The optimal convergence rate of the p-version of
the finite-element method. SIAM Journal on Numerical Analysis, vol. 24, no. 4, pages
750–776, 1987.
[Babuška & Suri 1987b] Ivo Babuška et Manil Suri. The h-p version of the finite element
method with quasiuniform meshes. RAIRO-Modélisation mathématique et analyse
numérique, vol. 21, no. 2, pages 199–238, 1987.
[Babuska & Szabo 1982] I. Babuska et B. Szabo. On the rates of convergence of the finite
element method. International Journal for Numerical Methods in Engineering, vol. 18,
no. 3, pages 323–341, 1982.
[Babuska et al. 1981] I. Babuska, B.A. Szabo et I.N. Katz. The p-version of the finite element
method. SIAM Journal on Numerical Analysis, vol. 18, no. 3, pages 515–545, 1981.
[Babuška et al. 1994] I. Babuška, T. Strouboulis, A. Marthur et C.S. Upadhyay. Pollution-error
in the h-version of the finite-element method and the local quality of a-posteriori error
estimators. Finite Elements in Analysis and Design, vol. 17, pages 273–321, 1994.
[Bailly et al. 1997] H. Bailly, D. Ménessier et C. Prunier. Le combustible nucléaire des réac-
teurs à eau sous pression et des réacteurs à neutrons rapides. Commissariat à l’Énergie
Atomique, 1997.
[Bank & Weiser 1985] R.E. Bank et A. Weiser. Some a posteriori error estimators for elliptic
partial differential equations. Mathematics of Computation, vol. 44, no. 170, pages
283–301, 1985.
[Bank & Welfert 1991] R.E. Bank et B.D. Welfert. A Posteriori Error Estimates for the Stokes
Problem. SIAM Journal on Numerical Analysis, vol. 28, no. 3, pages 591–623, 1991.
[Bank et al. 1983] R.E. Bank, A.H. Sherman et A. Weiser. Some refinement algorithms and
data structures for regular local mesh refinement. Scientific Computing, Applications of
Mathematics and Computing to the Physical Sciences, vol. 1, pages 3–17, 1983.
[Bao et al. 2012] G. Bao, G. Hu et D. Liu. An h-adaptive finite element solver for the calculations
of the electronic structures. Journal of Computational Physics, vol. 231, no. 14, pages
4967–4979, 2012.
[Barbié et al. 2014] L. Barbié, I. Ramière et F. Lebon. Strategies around the local defect correc-
tion multi-level refinement method for three-dimensional linear elastic problems. Com-
puters and Structures, vol. 130, pages 73–90, 2014.
Bibliographie 121

[Barbié et al. 2015] L. Barbié, I. Ramière et F. Lebon. An automatic multilevel refinement

technique based on nested local meshes for nonlinear mechanics. Computers & Structures,
vol. 147, pages 14–25, 2015.
[Barbié 2013] L. Barbié. Raffinement de maillage multi-grille local en vue de la simulation
3D du combustible nucléaire des Réacteurs à Eau sous Pression. Thèse, Aix-Marseille
Université, Octobre 2013.
[Barros et al. 2004] F.B. Barros, S.P.B. Proenca et C.S. de Barcellos. Generalized finite ele-
ment method in structural nonlinear analysis - a p-adaptive strategy. Computational
Mechanics, vol. 33, no. 2, pages 95–107, 2004.
[Beckett et al. 2001] G. Beckett, J.A. Mackenzie et M.L. Robertson. A moving mesh finite
element method for the solution of two-dimensional Stefan problems. Journal of Compu-
tational Physics, vol. 168, no. 2, pages 500–518, 2001.
[Bellenger & Coorevits 2005] E. Bellenger et P. Coorevits. Adaptive mesh refinement for the
control of cost and quality in finite element analysis. Finite Elements in Analysis and
Design, vol. 41, no. 15, pages 1413–1440, 2005.
[Belliard & Grandotto 2003] M. Belliard et M. Grandotto. Local zoom computation of two-
phase flows in steam generators using a local defect correction method. Numerical Heat
Transfer - Part A Applications, vol. 43, no. 2, pages 111–135, 2003.
[Belytschko & Tabbara 1993] T. Belytschko et M. Tabbara. H-adaptive finite-element methods
for dynamic problems, with emphasis on localization. International Journal for Numeri-
cal Methods in Engineering, vol. 36, no. 24, pages 4245–4265, 1993.
[Biotteau 2010] E. Biotteau. Stratégie multigrille et raffinement automatique à précision contro-
lée pour la dynamique transitoire non-linéaire. PhD thesis, Institut National des Sciences
Appliquées de Lyon, 2010.
[Boersma & Wriggers 1997] A Boersma et P Wriggers. An algebraic multigrid solver for finite
element computations in solid mechanics. Engineering Computations, vol. 14, no. 2,
pages 202–215, 1997.
[Bonnet & Frangi 2007] de Marc Bonnet et Attilio Frangi. Analyse des solides déformables par
la méthode des éléments finis. European Journal of Computational Mechanics/Revue
Européenne de Mécanique Numérique, vol. 16, no. 5, pages 667–668, 2007.
[Bouillard et al. 1996] P. Bouillard, J.F. Allard et G. Warzee. Superconvergent patch recovery
technique for the finite element method in acoustics. Communications in Numerical
Methods in Engineering, vol. 12, no. 9, pages 581–594, 1996.
[Boussetta et al. 2006] R. Boussetta, T. Coupez et L. Fourment. Adaptive remeshing based
on a posteriori error estimation for forging simulation. Computer Methods in Applied
Mechanics and Engineering, vol. 195, no. 48-49, pages 6626–6645, 2006.
[Brandt 1977] A. Brandt. Multi-level adaptive solutions to boundary-value problems. Mathe-
matics of Computation, vol. 31, pages 333–390, 1977.
[Brandt 1994] A. Brandt. Rigorous quantitative-analysis of Multigrid .1. Constant-coefficients
2-Level cycle with L2-norm. SIAM Journal on Numerical Analysis, vol. 31, no. 6, pages
1695–1730, 1994.
122 Bibliographie

[Cai & Li 2010] Z. Cai et R. Li. An h-adaptive mesh method for Boltzmann-BGK/hydrodynamics
coupling. Journal of Computational Physics, vol. 229, pages 1661–1680, 2010.
[Cai et al. 1995] Z. Cai, F. Le Gland et H. Zhang. An adaptive local grid refinement method for
nonlinear filtering. Rapport technique, INRIA, 1995.
[Cao et al. 2001] W. Cao, W. Huang et R.D. Russell. Comparison of two-dimensional r-adaptive
finite element methods using various error indicators. Mathematics and Computers in
Simulation, vol. 56, pages 127–143, 2001.
[Carstensen & Funken 1999] C. Carstensen et S. Funken. Fully Reliable Localized Error
Control in the FEM. SIAM Journal on Scientific Computing, vol. 21, no. 4, pages
1465–1484, 1999.
[CAST3M 2016] CAST3M. [Link], 2016.
[Cavin et al. 2005] P. Cavin, A. Gravouil, A.A. Lubrecht et A. Combescure. Efficient FEM
calculation with predefined precision through automatic grid refinement. Finite Elements
in Analysis and Design, vol. 41, pages 1043–1055, 2005.
[Chabrand et al. 1998] P. Chabrand, F. Dubois et M. Raous. Various numerical methods for sol-
ving unilateral contact problems with friction. Mathematical and Computer Modelling,
vol. 28, no. 4, pages 97–108, 1998.
[Ciarlet 1978] P.G. Ciarlet. The finite element method for elliptic problems, volume 4 of
Studies in Mathematics and its Applications. North-Holland Publishing Company, J.
Lions, G. Papanicolaou and R.T. Rockafellar édition, 1978.
[Clément 1975] Ph. Clément. Approximation by finite element functions using local regula-
rization. ESAIM : Mathematical Modelling and Numerical Analysis - Modélisation
Mathématique et Analyse Numérique, vol. 9, no. R2, pages 77–84, 1975.
[Cocu et al. 1996] M. Cocu, E. Pratt et M. Raous. Formulation and approximation of quasistatic
frictional contact. International Journal of Engineering Science, vol. 34, no. 7, pages
783–798, 1996.
[Cocu 1984] M. Cocu. Existence of solutions of Signorini problems with friction. International
Journal of Engineering Science, vol. 22, no. 5, pages 567–575, 1984.
[Code_Aster 2016] Code_Aster. [Link], 2016.
[Coorevits et al. 2000] P. Coorevits, P. Hild et J.P. Pelle. A posteriori error estimation for
unilateral contact with matching and non-matching meshes. Computer Methods in
Applied Mechanics and Engineering, vol. 186, no. 1, pages 65–83, 2000.
[Dari et al. 1999] E. Dari, R.G. Durán et C. Padra. Maximum norm error estimators for three-
dimensional elliptic problems. SIAM Journal on Numerical Analysis, vol. 37, no. 2,
pages 683–700, 1999.
[Dassault_Systèmes 2016] Dassault_Systèmes. 12.3.2 Error indicators, 2016.
[Delmas 2008] J. Delmas. Stratégies de contrôle d’erreur en calcul de structures industrielles.
PhD thesis, Université de Picardie Jules Verne, 2008.
[Delmas 2011] J. Delmas. Estimateurs d’erreur en résidu, 2011. Manuel de référence
Code_Aster.
Bibliographie 123

[Demkowicz et al. 1985] L. Demkowicz, P. Devloo et J.T. Oden. On a h-type mesh-refinement

strategy based on minimization of interpolation errors. Computer Methods in Applied
Mechanics and Engineering, vol. 53, no. 1, pages 67–89, 1985.
[Demkowicz et al. 1989] L. Demkowicz, J.T. Oden, W. Rachowicz et O. Hardy. Toward a
universal h-p adaptive finite-element strategy - 1. Constrained approximation and data
structure. Computer Methods in Applied Mechanics and Engineering, vol. 77, no. 1-2,
pages 79–112, 1989.
[Düster & Rank 2001] A. Düster et E. Rank. The p-version of the finite element method compa-
red to an adaptive h-version for the deformation theory of plasticity. Computer Methods
in Applied Mechanics and Engineering, vol. 190, pages 1925–1935, 2001.
[Düster et al. 2007] A. Düster, A. Niggl et E. Rank. Applying the hp-d version of the FEM to
locally enhance dimensionally reduced models. Computer Methods in Applied Mechanics
and Engineering, vol. 196, pages 3524–3533, 2007.
[Duvaut & Lions 1976] G. Duvaut et J.L. Lions. Inequalities in mechanics and physics, volume
219. Springer-Verlag, 1976.
[Edwards et al. 1993] M.G. Edwards, J.T. Oden et L. Demkowicz. An hr-adaptive approximate
Riemann solver for the Euler equations in 2 dimensions. SIAM Journal on Scientific
Computing, vol. 14, no. 1, pages 185–217, 1993.
[Ehlers et al. 2002] W. Ehlers, M. Ammann et S. Diebels. h-Adaptive FE methods applied
to single-and multiphase problems. International Journal for Numerical Methods in
Engineering, vol. 54, no. 2, pages 219–239, 2002.
[Fedorenko 1961] R.P. Fedorenko. A relaxation method for solving elliptic difference equa-
tions. USSR Computational Mathematics and Mathematical Physics, pages 1092–1096,
1961.
[Ferket & Reusken 1996a] P.J.J. Ferket et A. Reusken. A Finite Difference Discretization Me-
thod for Elliptic Problems Composite Grids. Computing, vol. 56, no. 4, pages 343–370,
1996.
[Ferket & Reusken 1996b] P.J.J. Ferket et A.A. Reusken. Further analysis of the local defect
correction method. Computing, vol. 56, no. 2, pages 117–139, 1996.
[Ferket 1996] P.J.J. Ferket. The finite difference based fast adaptive composite grid method.
Numerical Linear Algebra with Applications, vol. 3, no. 5, pages 391–411, 1996.
[Fichera 1973] Gaetano Fichera. Boundary value problems of elasticity with unilateral
constraints. In Linear Theories of Elasticity and Thermoelasticity, pages 391–424.
Springer, 1973.
[Fish & Guttal 1996] J. Fish et R. Guttal. The s-version of finite element method for laminated
composites. International Journal for Numerical Methods in Engineering, vol. 39,
no. 21, pages 3641–3662, 1996.
[Fish & Markolefas 1993] J. Fish et S. Markolefas. Adaptive s-method for linear elastostatics.
Computer Methods in Applied Mechanics and Engineering, vol. 104, pages 363–396,
1993.
124 Bibliographie

[Fish & Markolefas 1994] J. Fish et S. Markolefas. Adaptive global-local refinement strategy
based on the interior error-estimatesd of the h-method. International Journal for Nume-
rical Methods in Engineering, vol. 37, no. 5, pages 827–838, 1994.
[Fish 1992] J. Fish. The s-version of the finite-element method. Computers and Structures,
vol. 43, no. 3, pages 539–547, 1992.
[Gallimard 2009] L. Gallimard. A constitutive relation error estimator based on traction-free
recovery of the equilibrated stress. International Journal for Numerical Methods in
Engineering, vol. 78, no. 4, pages 460–482, 2009.
[Garlick 1973] A. Garlick. Fracture of Zircaloy cladding under simulated power ramp conditions.
Journal of Nuclear Materials, vol. 49, pages 209–224, 1973.
[Ghosh & Manna 1993] S. Ghosh et S.K. Manna. r-adapted arbitrary Lagrangian-Eulerian
finite-element method in metal-forming simulation. Journal of Materials Engineering
and Performance, vol. 2, no. 2, pages 271–282, 1993.
[Glowinski et al. 1976] R. Glowinski, J.-L. Lions et R. Tremolieres. Analyse numerique des
inequations variationnelles. Dunod, 1976.
[Gratsch & Bathe 2005] T. Gratsch et K.J. Bathe. Review - A posteriori error estimation tech-
niques in practical finite element analysis. Computers and Structures, vol. 83, pages
235–265, 2005.
[Grego 1995] Laurence Grego. Méthodes multiniveaux pour des problèmes de contact unilatéral
avec frottement. PhD thesis, Université Aix-Marseille 1, 1995.
[Guerin et al. 2013] Y. Guerin, D. Parrat et J. Noirot. Nuclear Fuels for Light Water Reactors
and Fast Reactors. International School In Nuclear Engineering, Novembre 2013.
[Guo & Babuška 1986] B Guo et I Babuška. The hp version of the finite element method.
Computational Mechanics, vol. 1, no. 1, pages 21–41, 1986.
[Hackbusch 1984] W. Hackbusch. Local Defect Correction Method and Domain Decomposition
Techniques. Computing Suppl. Springer-Verlag, vol. 5, pages 89–113, 1984.
[Hackbusch 1985] W. Hackbusch. Multi-Grid Methods and Applications. Numéro 4 de Sprin-
ger Series in Computational Mathematics. Springer-Verlag, 1985.
[Han & Wriggers 2000] C.S. Han et P. Wriggers. An h-adaptive method for elasto-plastic shell
problems. Computer Methods in Applied Mechanics and Engineering, vol. 189, pages
651–671, 2000.
[Hart et al. 1986] L. Hart, S. McCormick, A. O’gallagher et J. Thomas. The fast adaptive
composite-grid method (FAC) : algorithms for advanced computers. Applied Mathematics
and Computation, vol. 19, no. 1, pages 103–125, 1986.
[Heegaard & Curnier 1993] J.H. Heegaard et A. Curnier. An augmented Lagrangian method
for discrete large-slip contact problems. International Journal for Numerical Methods
in Engineering, vol. 36, no. 4, pages 569–593, 1993.
[Huang & Russell 2010] W. Huang et R.D. Russell. Adaptive moving mesh methods, volume
174. Springer Science & Business Media, 2010.
[Jernkvist 1995] L. O. Jernkvist. A model for predicting pellet-cladding interaction-induced fuel
rod failure. Nuclear Engineering and Design, vol. 156, pages 393–399, 1995.
Bibliographie 125

[Kelly et al. 1983] D.W. Kelly, J.P. De S.R. Gago, O.C. Zienkiewicz et I. Babuška. A posteriori
error analysis and adaptive processes in the finite element method : Part I—error analysis.
International Journal for Numerical Methods in Engineering, vol. 19, no. 11, pages
1593–1619, 1983.
[Khadra et al. 1996] K. Khadra, Ph. Angot, J.P. Caltagirone et P. Morel. Concept de zoom
adaptatif en architecture multigrille locale ; étude comparative des méthodes L.D.C., F.A.C.
et F.I.C. RAIRO - Modélisation Mathématique et Analyse Numérique, vol. 30, no. 1,
pages 39–82, 1996.
[Kikuchi & Oden 1988] N. Kikuchi et J.T. Oden. Contact problems in elasticity : a study of
variational inequalities and finite element methods, volume 8. SIAM, 1988.
[Kikuchi & Song 1981] N. Kikuchi et Y. J. Song. Penalty/finite-element approximations of a
class of unilateral problems in linear elasticity. Quarterly of Applied Mathematics, pages
1–22, 1981.
[Klarbring & Björkman 1988] A. Klarbring et G. Björkman. A mathematical programming
approach to contact problems with friction and varying contact surface. Computers &
Structures, vol. 30, no. 5, pages 1185–1198, 1988.
[Klarbring 1986] A. Klarbring. A mathematical programming approach to three-dimensional
contact problems with friction. Computer Methods in Applied Mechanics and Enginee-
ring, vol. 58, no. 2, pages 175–200, 1986.
[Kramer et al. 2009] W. Kramer, H.J.H. Clercx, R.M.M. Mattheij et R. Minero. A finite volume
local defect correction method for solving the transport equation. Computers and Fluids,
vol. 38, no. 3, pages 533–543, 2009.
[Kubatko et al. 2009] E.J. Kubatko, S. Bunya, C. Dawson et J.J. Westerink. Dynamic p-
adaptive Runge-Kutta discontinuous Galerkin methods for the shallow water equations.
Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 21-26, pages
1766–1774, 2009.
[Kuss & Lebon 2011] F. Kuss et F. Lebon. Error estimation and mesh adaptation for Signorini-
Coulomb problems using E-FEM. Computers and Structures, vol. 89, no. 11-12, pages
1148–1154, 2011.
[Ladevèze & Leguillon 1983] P. Ladevèze et D. Leguillon. Error estimate procedure in the
finite-element method and applications. SIAM Journal on Numerical Analysis, vol. 20,
pages 485–509, 1983.
[Ladevèze & Pelle 2005] P. Ladevèze et J.P. Pelle. Mastering calculations in linear and nonli-
near mechanics. Springer, 2005.
[Ladevèze et al. 1986] P. Ladevèze, G. Coffignal et J.P. Pelle. Accuracy of elastoplastic and
dynamic analysis. In I. Babuska, O.C. Zienkiewicz, J.P. Gago et Oliveira E.R., edi-
teurs, Accuracy estimates and adaptive refinements in Finite Element computations,
chapitre 11, pages 181–203. Wiley, New York, 1986.
[Ladevèze et al. 1991] P. Ladevèze, J.P. Pelle et Ph. Rougeot. Error estimation and mesh
optimization for classical finite elements. Engineering Computations, vol. 8, no. 1,
pages 69–80, 1991.
126 Bibliographie

[Ladevèze 1975] P. Ladevèze. Comparaison de modèles de milieux continus. PhD thesis,

Université Pierre et Marie Curie, Paris, 1975.
[Ladevèze et al. 2010] P. Ladevèze, L. Chamoin et É. Florentin. A new non-intrusive technique
for the construction of admissible stress fields in model verification. Computer Methods
in Applied Mechanics and Engineering, vol. 199, no. 9–12, pages 766 – 777, 2010.
[Lax & Milgram 1954] P.D. Lax et A.N. Milgram. IX. Parabolic equations. Contributions to
the Theory of Partial Differential Equations, no. 33, page 167, 1954.
[Lebon & Raous 1992] F. Lebon et M. Raous. Multibody contact problem including friction in
structure assembly. Computers & Structures, vol. 43, no. 5, pages 925 – 934, 1992.
[Lebon et al. 2007] F. Lebon, M. Raous et I. Rosu. Multigrid Methods for Unilateral Contact
Problems with Friction. In P. Wriggers et U. Nackenhorst, editeurs, IUTAM Symposium
on Computational Contact Mechanics, pages 1–16. Springer, 2007. 5-8 November,
Hannover, Germany.
[Lebon 1995] F. Lebon. Two-grid method for regularized frictional elastostatics problems. En-
gineering Computations, vol. 12, pages 657–664, 1995.
[Lebrun-Grandié et al. 2011] D. Lebrun-Grandié, J. Ragusa et B. Turcksin. Adaptive multimesh
hp-FEM for a coupled neutronics and nonlinear heat conduction problem. In International
Conference on Mathematics and Computational Methods Applied to Nuclear Science
and Engineering. American Nuclear Society, 2011. 8-12 May, Rio de Janeiro, Brazil.
[Lemke 1980] C.E. Lemke. A survey of complementarity theory. In Variational Inequalities
and Complementarity Problems, pages 213–235. John Wiley, 1980.
[Li et al. 1995] L. Li, P. Bettess, J.W. Bull, T. Bond et I. Applegarth. Theoretical formulations
for adaptive finite element computations. Communications in Numerical Methods in
Engineering, vol. 11, no. 10, pages 857–868, 1995.
[Lubrecht et al. 1987] A.A. Lubrecht, W.E. Ten Napel et R. Bosma. Multigrid, an alternative
method of solution for two-dimensional elastohydrodynamically lubricated point contact
calculations. Journal of Tribology, vol. 109, no. 3, pages 437–443, 1987.
[Machiels et al. 2000] L. Machiels, Y. Maday et A.T. Patera. A “flux-free” nodal Neumann sub-
problem approach to output bounds for partial differential equations. Comptes Rendus
de l’Académie des Sciences - Series I - Mathematics, vol. 330, no. 3, pages 249 – 254,
2000.
[Marchal et al. 2009] N. Marchal, C. Campos et C. Garnier. Finite element simulation of
Pellet-Cladding Interaction (PCI) in nuclear fuel rods. Computational Materials Science,
vol. 45, pages 821–826, 2009.
[Marelle & Thouvenin 2011] V. Marelle et G. Thouvenin. Alcyone V1.2.3. : Notice d’utilisation.
Rapport technique NT 11-012, CEA/DEN/CAD/DEC/SESC/LSC, 2011.
[McCormick & Thomas 1986] S.F. McCormick et J. Thomas. The fast adaptive composite grid
(FAC) method for elliptic-equations. Mathematics of Computation, vol. 46, no. 174,
pages 439–456, 1986.
[McCormick 1984] S.F. McCormick. Fast Adaptive Composite Grid (F.A.C.) Methods : theory
for the variational case. Computing Suppl. Springer-Verlag, vol. 5, pages 115–121,
1984.
Bibliographie 127

[Michel et al. 2008a] B. Michel, J. Sercombe et G. Thouvenin. A new phenomenological crite-

rion for pellet-cladding interaction rupture. Nuclear Engineering and Design, vol. 238,
pages 1612–1628, 2008.
[Michel et al. 2008b] B. Michel, J. Sercombe, G. Thouvenin et R. Chatelet. 3D fuel cracking
modelling in pellet cladding mechanical interaction. Engineering Fracture Mechanics,
vol. 75, pages 3581–3598, 2008.
[Michel et al. 2013] B. Michel, C. Nonon, J. Sercombe, F. Michel et V. Marelle. Simulation
of pellet-cladding interaction with the PLEIADES fuel performance software environment.
Nuclear Technology, vol. 182, 2013.
[Minero et al. 2006] R. Minero, M.J.H. Anthonissen et R.M.M. Mattheij. A local defect correc-
tion technique for time-dependent problems. Numerical Methods for Partial Differential
Equations, vol. 22, pages 128–144, 2006.
[Mitchell 1991] W.F. Mitchell. Adaptive refinement for arbitrary finite-element spaces with
hierarchical bases. Journal of Computational and Applied Mathematics, vol. 36, no. 1,
pages 65–78, 1991.
[Mougel et al. 2004] C. Mougel, B. Verhaeghe, C. Verdeau, S. Lansiart, S. Béguin et B. Julien.
Power ramping in the OSIRIS reactor : data base analysis for standard UO2 fuel with
Zy-4 cladding. In International Seminar on Pellet-Clad Interaction on Water Reactor
Fuels, March 2004. 9-11 March, Aix-en-Provence, France.
[Nambiar & Lawrence 1992] R.V. Nambiar et K.L. Lawrence. The Zienkiewicz–Zhu error esti-
mator for multiple material problems. Communications in Applied Numerical Methods,
vol. 8, no. 4, pages 273–277, 1992.
[Nefedov & Mattheij 2002] V. Nefedov et R.M.M. Mattheij. Local Defect Correction with Dif-
ferent Grid Types. Numerical Methods for Partial Differential Equations, vol. 18, no. 4,
pages 545–468, 2002.
[Nochetto 1995] R.H. Nochetto. Pointwise a posteriori error estimates for elliptic problems
on highly graded meshes. Mathematics of Computation, vol. 64, no. 209, pages 1–22,
1995.
[Nonon et al. 2004] C. Nonon, S. Lansiart, C. Struzik, D. Plancq, S. Martin, G.M. Decroix,
O. Rambouille, S. Beguin et B. Julien. Differential PCI behaviour of PWR fuel rods
under transient conditions. In International Topical Meeting on LWR Fuel Performance.
American Nuclear Society, 2004. 19-22 September, Orlando, Florida.
[Oden & Demkowicz 1991] J.T. Oden et L. Demkowicz. H-p adaptive finite-element methods
in computational fluid-dynamics. Computer Methods in Applied Mechanics and Engi-
neering, vol. 89, no. 1-3, pages 11–40, 1991.
[Oden et al. 1994] J.T. Oden, W. Wu et M. Ainsworth. An a posteriori error estimate for finite
element approximations of the Navier-Stokes equations. Computer Methods in Applied
Mechanics and Engineering, vol. 111, no. 1–2, pages 185 – 202, 1994.
[Oñate & Bugeda 1993] E. Oñate et G. Bugeda. A study of mesh optimality criteria in adaptive
finite element analysis. Engineering computations, vol. 10, no. 4, pages 307–321, 1993.
[Panagiotopoulos 1985] P.D. Panagiotopoulos. Inequality problems in mechanics and appli-
cations : Convex and nonconvex energy functions. Birkhaüser, Basel, 1985.
128 Bibliographie

[Pares et al. 2006] N. Pares, P. Díez et A. Huerta. Subdomain-based flux-free a posteriori error
estimators. Computer Methods in Applied Mechanics and Engineering, vol. 195, pages
297–323, 2006.
[Pled et al. 2011] F. Pled, L. Chamoin et P. Ladevèze. On the techniques for constructing
admissible stress fields in model verification : Performances on engineering examples.
International Journal for Numerical Methods in Engineering, vol. 88, no. 5, pages
409–441, 2011.
[Prudhomme et al. 2004] S. Prudhomme, F. Nobile, L. Chamoin et J.T. Oden. Analysis of a
subdomain-based error estimator for finite element approximations of elliptic problems.
Numerical Methods for Partial Differential Equations, vol. 20, no. 2, pages 165–192,
2004.
[Rachowicz & Zdunek 2005] W. Rachowicz et A. Zdunek. An hp-adaptive finite element me-
thod for scattering problems in computational electromagnetics. International Journal
for Numerical Methods in Engineering, vol. 62, no. 9, pages 1226–1249, 2005.
[Rachowicz et al. 1989] W. Rachowicz, J.T. Oden et L. Demkowicz. Toward a universal h-p
adaptive finite-element strategy - 3. Design of h-p meshes. Computer Methods in Applied
Mechanics and Engineering, vol. 77, no. 1-2, pages 181–212, 1989.
[Rachowicz 1997] W. Rachowicz. An anisotropic h-adaptive finite element method for compres-
sible Navier-Stokes equations. Computer Methods in Applied Mechanics and Enginee-
ring, vol. 146, pages 231–252, 1997.
[Ramière 2008] I. Ramière. Convergence analysis of the Q1 -finite element method for ellip-
tic problems with non-boundary-fitted meshes. International Journal for Numerical
Methods in Engineering, vol. 75, no. 9, pages 1007–1052, 2008.
[Ramm et al. 2003] E. Ramm, E. Rank, R. Rannacher, K. Schweizerhof, E. Stein, W. Wendland,
G. Wittum, P. Wriggers et W. Wunderlich. Error-controlled adaptive finite elements in
solid mechanics. John Wiley & Sons, 2003.
[Raous et al. 1988] M. Raous, P. Chabrand et F. Lebon. Numerical methods for frictional
contact problems and applications. Journal of Theoretical and Applied Mechanics,
vol. 7, pages 111–128, 1988.
[Raous 1999] M. Raous. Quasistatic Signorini problem with Coulomb friction and coupling to
adhesion. In New Developments in Contact Problems, pages 101–178. Springer, 1999.
[Raviart & Thomas 1998] P-A. Raviart et J-M. Thomas. Introduction à l’analyse numérique
des équations aux dérivées partielles. P.G. Ciarlet and J.L. Lions, Dunod édition, 1998.
[Ritter et al. 2010] D. Ritter, M. Sturmer et U. Rude. A fast-adaptive composite grid algorithm
for solving the free-space Poisson problem on the cell broadband engine. Numerical
Linear Algebra with Applications, vol. 17, no. 2-3, pages 291–305, 2010.
[Rivara 1984] M. Rivara. Mesh refinement processes based on the generalized bisection of
simplices. SIAM Journal on Numerical Analysis, vol. 21, no. 3, pages 604–613, 1984.
[Roberts 1978] G. Roberts. The concentration of stress in cladding produced by the expansion
of cracked fuel pellets. Nuclear Engineering and Design, vol. 47, pages 257–266, 1978.
Bibliographie 129

[Sercombe et al. 2012] J. Sercombe, I. Aubrun et C. Nonon. Power ramped cladding stresses
and strains in 3D simulations with burnup-dependent pellet-clad friction. Nuclear Engi-
neering and Design, vol. 242, pages 164–181, 2012.
[Sercombe et al. 2013] J. Sercombe, J. Julien, F. Michel, B. Michel et E. Fédérici. 3D modelling
of strain concentration due to PCI within the fuel code ALCYONE. In Proceedings of Top
Fuel 2013. American Nuclear Society, September 15-19 2013. Charlotte, NC.
[Simo & Laursen 1992] J.C. Simo et T.A. Laursen. An augmented Lagrangian treatment of
contact problems involving friction. Computers & Structures, vol. 42, no. 1, pages
97–116, 1992.
[Solin et al. 2010] P. Solin, J. Ceverny, L. Dubcova et D. Andrs. Monolithic Discretization of
Linear Thermoelasticity Problems via Adaptative Multimesh hp-FEM. Journal of Compu-
tational and Applied Mathematics, vol. 234, pages 2350–2357, 2010.
[Strouboulis & Haque 1992] T. Strouboulis et K.A. Haque. Recent experiences with error esti-
mation and adaptivity - 2. Error estimation for h-adaptive approximations on grids of
triangles and quadrilaterals. Computer Methods in Applied Mechanics and Enginee-
ring, vol. 100, no. 3, pages 359–430, 1992.
[Sun & Zamani 1992] W. Sun et N.G. Zamani. An adaptive h-r boundary element algorithm
for the Laplace equation. International Journal for Numerical Methods in Engineering,
vol. 33, no. 3, pages 537–552, 1992.
[Thomas et al. 1986] J.W. Thomas, R. Schweitzer, M. Heroux, S. McCormick et A.M. Thomas.
Application of the fast adaptive composite grid method to computational fluid dynamics.
In Tenth International Conference on Numerical Methods in Fluid Dynamics, pages
606–611. Springer, 1986.
[Verfürth 1996] R. Verfürth. A review of a posteriori error estimation and adaptive mesh-
refinement techniques, volume 1. Wiley-Teubner Chichester, 1996.
[Verfürth 1999] R. Verfürth. A review of a posteriori error estimation techniques for elasticity
problems. Computer Methods in Applied Mechanics and Engineering, vol. 176, no. 1,
pages 419–440, 1999.
[Wappler 1999] J. U. Wappler. Die lokale Defektkorrekturmethode zur adaptiven Diskretisie-
rung elliptischer Differentialgleichungen mit finiten Elementen. PhD thesis, Christian-
Albrechts-Universität, Kiel, 1999. in German.
[Watremetz et al. 2007] B. Watremetz, M.C. Baietto-Dubourg et A.A. Lubrecht. 2D thermo-
mechanical contact simulations in a functionally graded material : A multigrid-based
approach. Tribology International, vol. 40, pages 754–762, 2007.
[Wriggers & Boersma 1998] P. Wriggers et A. Boersma. A parallel algebraic multigrid solver
for problems in solid mechanics discretisized by finite elements. Computers & Structures,
vol. 69, no. 1, pages 129–137, 1998.
[Wriggers & Scherf 1998] P. Wriggers et O. Scherf. Different a posteriori error estimators
and indicators for contact problems. Mathematical and Computer Modelling, vol. 28,
no. 4-8, pages 437–447, 1998.
130 Bibliographie

[Wriggers & Zavarise 1993] P. Wriggers et G. Zavarise. Application of augmented Lagrangian

techniques for non-linear constitutive laws in contact interfaces. Communications in
Numerical Methods in Engineering, vol. 9, no. 10, pages 815–824, 1993.
[Wriggers 2006] P. Wriggers. Computational contact mechanics. Springer Berlin Heidelberg,
2006.
[Yue & Robbins 2005] Z. Yue et D.H. Robbins. Adaptive superposition of finite element meshes
in elastodynamic problems. International Journal for Numerical Methods in Enginee-
ring, vol. 63, no. 11, pages 1604–1635, 2005.
[Yue & Robbins 2007] Z. Yue et D.H. Robbins. Adaptive superposition of finite element meshes
in non-linear transient solid mechanics problems. International Journal for Numerical
Methods in Engineering, vol. 72, pages 1063–1094, 2007.
[Zhu 1997] J. Zhu. A posteriori error estimation—the relationship between different procedures.
Computer Methods in Applied Mechanics and Engineering, vol. 150, no. 1, pages
411–422, 1997.
[Zienkiewicz & Zhu 1987] O.C. Zienkiewicz et J.Z. Zhu. A simple error estimator and adap-
tive procedure for practical engineering analysis. International Journal for Numerical
Methods in Engineering, vol. 24, pages 337–357, 1987.
[Zienkiewicz & Zhu 1992a] O.C. Zienkiewicz et J.Z. Zhu. The superconvergent patch recovery
and a posteriori error estimation. Part I : The recovery technique. International Journal
for Numerical Methods in Engineering, vol. 33, pages 1331–1364, 1992.
[Zienkiewicz & Zhu 1992b] O.C. Zienkiewicz et J.Z. Zhu. The superconvergent patch recovery
and a posteriori error estimation. Part II : Error estimates and adaptivity. International
Journal for Numerical Methods in Engineering, vol. 33, pages 1365–1382, 1992.
[Zienkiewicz et al. 1983] O.C. Zienkiewicz, J.D.S. Gago et D.W. Kelly. The hierarchical concept
in finite element analysis. Computers & Structures, vol. 16, no. 1-4, pages 53–65, 1983.
[Zienkiewicz et al. 2005] O.C. Zienkiewicz, J. Zhu et R.L. Taylor. The finite element method :
its basis and fundamentals. Butterworth-Heinemann, 2005.