# Plasticité des structures

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Plasticité des structures. École d'ingénieur. Ecole Centrale de Nantes, France. 2015. cel-02059317

HAL Id: cel-02059317

https://cel.hal.science/cel-02059317v1

Submitted on 6 Mar 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Thomas Heuzé

# Plasticité des structures

<figure><img src="image_2.png" /></figure>

Centrale Nantes
École Centrale de Nantes — 2015

Ce document est sous licence Creative Commons BY-NC-SA 4.0 France :

— attribution ;

— pas d'utilisation commerciale ;

— partage dans les mêmes conditions.

http://creativecommons.org/licenses/by-nc-sa/4.0/deed.fr

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>
# Table des matières

## Chapitre 1 Élastoplasticité en petites transformations

### 1.1 Phénoménologie globale 1

1.1.1 Comportement élastoplastique uniaxial : l'essai de traction 1

1.1.2 Sollicitation multiaxiale 3

### 1.2 Notions sur les origines physiques 4

### 1.3 Critères de plasticité 6

1.3.1 Critère de von Mises 8

1.3.2 Critère de Tresca 9

1.3.3 Comparaison des critères de von Mises et Tresca 10

### 1.4 Modélisation de l'écrouissage 11

1.4.1 Écrouissage isotrope 12

1.4.2 Écrouissage cinématique 13

### 1.5 Règle d'écoulement plastique 14

1.5.1 Partition de la déformation 14

1.5.2 Principe de Hill 15

1.5.3 Calcul du multiplicateur plastique 17

1.5.4 Postulat de Drucker-Ilyushin et principe de Hill 18

## Chapitre 2 Thermodynamique des milieux continus

### 2.1 Rappels de thermodynamique 21

2.1.1 Premier principe de la thermodynamique 21

2.1.2 Second principe de la thermodynamique 22

2.1.3 Écriture des lois de comportement 24

### 2.2 Méthode de l'état local 24

2.2.1 Potentiel thermodynamique et lois d'état 25

2.2.2 Potentiel de dissipation et lois complémentaires 26

2.2.3 Thermodynamique des processus irréversibles 27

### 2.3 Plasticité associée 28

2.3.1 Principe de la maximisation de la dissipation intrinsèque 28

2.3.2 Particularisation au critère de von Mises 28

2.3.3 Matériaux élastique et parfaitement plastique 29

2.3.4 Lois à écrouissage isotrope 30

2.3.5 Lois à écrouissage cinématique linéaire 31

### 2.4 Plasticité non associée 32

2.4.1 Calcul du multiplicateur plastique 33

2.4.2 Particularisation 33
<figure><img src="image_8.png" /></figure>

## Chapitre 3 Calcul des structures élastoplastiques

3.1 Un problème d'évolution 37

3.2 Formulation forte du problème 38

3.3 Résultats sur l'évolution des systèmes élastoplastiques 39

3.4 Problèmes d'élastoplasticité en petites transformations 40

3.4.1 Structure treillis à trois barres 40

3.4.2 Exercices sur les barres en élastoplasticité 44

3.4.3 Plasticité dans les poutres en flexion simple 45

3.5 État résiduel d'un système élastoplastique 50

3.5.1 Incompatibilité des déformations plastiques et contraintes résiduelles 50

3.5.2 Paramètres généralisés de chargement 52

3.5.3 Comportement irréversible du système 52

3.5.4 Écrouissage de structure en variables généralisées 54

## Chapitre 4 Simulation numérique des structures élastoplastiques

4.1 Simulation numérique en calcul des structures 55

4.2 Méthode des éléments finis en mécanique non-linéaire 55

4.2.1 Formulation faible du problème élastoplastique 55

4.2.2 Discrétisation éléments finis et équations discrètes 56

4.2.3 Résolution du problème stationnaire 57

4.3 Intégration des lois de comportement élastoplastiques 60

4.3.1 Modèle de plasticité de von Mises avec écrouissage isotrope linéaire 60

4.3.2 Algorithmes de type Return Mapping 60

## Annexe A Travaux dirigés

A.1 Plasticité monodimensionnelle 71

A.1.1 Modèle de Saint-Venant 71

A.1.2 Plasticité en chargement monotone et cyclique 71

A.1.3 Modèle de Saint-Venant généralisé 72

A.2 Loi de comportement élastoplastique 72

A.3 Torsion élastoplastique d'un arbre cylindrique plein 73

A.4 Implémentation de lois de comportement à seuil 74

A.4.1 Intégration numérique temporelle 74

A.4.2 Implémentation numérique d'un modèle de plasticité 74

## Annexe B Devoirs surveillés

B.1 Sphère sous pression en élastoplasticité [2012-2013] 77

B.2 Analyse d'un système de barres en traction en élastoplasticité [2013-2014] 78

B.3 Implémentation d'une loi de comportement élastoplastique [2014-2015] 80

## Annexe C Projet numérique

C.1 Système unidimensionnel 83

C.1.1 Développement du code de calcul par éléments finis 83
<figure><img src="image_9.png" /></figure>

<table><tr><td>C.1.2</td><td>Application : pâle de turbine soumise à un chargement centrifuge</td><td>87</td></tr><tr><td>C.2</td><td>Système bidimensionnel axisymétrique</td><td>88</td></tr><tr><td>C.2.1</td><td>Développement du code de calcul par éléments finis</td><td>88</td></tr><tr><td>C.2.2</td><td>Application : sphère sous pression interne</td><td>92</td></tr><tr><td>C.3</td><td>Composantes de contrainte en coordonnées sphériques</td><td>95</td></tr></table>

Bibliographie
# Introduction

Les exigences de compétitivité rendant très aigus les problèmes d’économie de matière, il est bien souvent nécessaire d’optimiser le choix des matériaux pour le dimensionnement des structures. Aussi pour une application donnée, il importe de caractériser et de modéliser le comportement du matériau utilisé dans ses conditions d’emploi. Ces exigences conduisent à devoir effectuer des analyses plus fines et à utiliser des modèles de comportement plus complexes, lesquels intègrent généralement des effets non linéaires.

La théorie de la plasticité traite des déformations irréversibles indépendantes du temps, elle est basée sur des mécanismes physiques intervenant dans les métaux et alliages mettant en jeu des mouvements de dislocations dans un réseau cristallin sans influence de phénomènes visqueux ni présence de décohésion endommageant la matière. Une des caractéristiques de la plasticité est qu’elle n’apparaît qu’une fois un seuil de charge atteint. Associée à cet effet de seuil, matérialisé par le domaine d’élasticité, on observe alors expérimentalement l’irréversibilité d’une partie de la déformation générée, notamment en faisant une décharge. Aussi, le comportement de l’élément de matière dépend du trajet de chargement subi, et donc des effets d’histoire.

Ce type de modèle de comportement utilisé en calcul des structures conduit à des analyses non-linéaires qui sont bien sûr plus complexes et plus coûteuses que les calculs menés en élasticité linéaire, cependant elles permettent de mieux s’approcher du comportement réel du système considéré et donc d’effectuer un dimensionnement plus strict des structures, notamment en vue de l’adapter aux sollicitations qu’elle subira (choix du matériau, optimisation de la forme de la pièce, etc) ou d’améliorer les prévisions de durée de vie de la structure.

L’objet de ce cours est d’étudier la « plasticité dans les structures », et donc de faire une présentation du comportement élastoplastique des solides et d’aller jusqu’à sa prise en compte dans un calcul de structure, et à l’implémentation de modèles élastoplastiques dans un code éléments finis. Le cours est organisé comme suit :

- le premier chapitre est consacré à l’étude du comportement élastoplastique des solides en petites transformations. Après une présentation phénoménologique de la plasticité et quelques notions sur ses origines à l’échelle microscopique, les équations de base de la plasticité (critère, loi d’écoulement) sont présentées;

- le second chapitre a pour objet de formaliser les équations de la plasticité, présentées au premier chapitre sur la base de principes, dans un cadre thermodynamique. On montre notamment que ce cadre s’avère puissant pour la construction de lois de comportement, et conduit notamment au modèle des matériaux standards généralisés (dans lequel la plasticité s’insère);

- le troisième chapitre traite du calcul des structures élastoplastiques. On y présente la formulation du problème aux limites avec un comportement élastoplastique, ainsi que quelques exemples de solutions analytiques de problèmes élastoplastiques en petites transformations;

- enfin, le quatrième chapitre s’intéresse à la simulation numérique des structures élastoplastiques. Après un bref rappel de la stratégie de résolution généralement suivi pour la résolution d’un problème non linéaire avec la méthode des éléments finis, l’attention est portée sur les algorithmes d’intégration d’une relation de comportement élastoplastique.
<figure><img src="image_10.png" /></figure>
# Élastoplasticité en petites transformations

D'une façon générale, les matériaux dans lesquels se développe la plasticité sont les matériaux métalliques caractérisés, à température ambiante, par une bonne ductilité, c'est-à-dire une aptitude à se déformer beaucoup avant rupture. Les déformations irréversibles générées par plasticité peuvent être exploitées, par exemple dans la mise en forme des métaux (laminage, forgeage, emboutissage...), ce type de procédé représentant en effet un mode de fabrication majeur des pièces mécaniques industrielles, ou être indésirables, auquel cas la structure en service reste saine en se déformant uniquement dans son domaine élastique.

Ce chapitre détaille le modèle de comportement élastoplastique des solides restreint au cadre des petites transformations. Ce modèle, comme bien d'autres en mécanique des solides, a été élaboré à partir de constations expérimentales, dont nous rappellerons les enseignements.

Comme tout modèle, le solide élastoplastique classique présente un certain domaine de validité dont il est fort utile de tracer le contour. Par exemple, la plasticité laisse de côté d'autres phénomènes physiques comme le vieillissement du matériau ou encore les effets de viscosité. En conséquence, la réponse du matériau à une sollicitation donnée est indépendante du temps physique. On verra que, bien que des dérivées temporelles interviennent dans les équations décrivant les lois d'évolution de certaines variables du modèle élastoplastique, ce temps n'a rien de physique (on parle de temps cinématique) et permet juste de décrire le séquencement des événements et donc les effets d'histoire de chargement subis par l'élément de matière. Qui plus est, les transformations considérées sont lentes, et les modèles développés valides dans le cas de chargements quasi-statiques et sans effet de la température.

## 1.1 Phénoménologie globale

Le comportement élastoplastique des solides est mis en évidence au travers d'essais de caractérisation mécaniques. Ces essais sont réalisés sur des éprouvettes dont la géométrie permet de remonter au champ de contrainte dans au moins une zone de cette éprouvette, qu'on apporte la zone utile. L'essai le plus simple à mettre en œuvre, et le plus utilisé, est l'essai de traction. D'autres essais génèrent des chargements multiaxiaux et auront pour but de qualifier le comportement multiaxial du matériau.

### 1.1.1 Comportement élastoplastique uniaxial : l'essai de traction

L'essai de traction consiste à tirer sur une barre, et permet d'extraire la courbe contrainte-déformation unidimensionnelle du matériau illustrée sur la figure 1.1. On observe au cours d'un chargement monotone croissant que la courbe présente deux phases :

**Première phase** Elle est définie par $\sigma < \sigma_0$. La relation entre $\sigma$ et $\varepsilon$ est linéaire. Si le chargement n'a pas dépassé la limite de contrainte $\sigma_0$, la décharge repassera par le même chemin :
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 1.1 – Essai de traction</figcaption></figure>

il y a réversibilité ; c'est la définition de l'élasticité. La relation de comportement s'écrit alors :

$$
\sigma = E\epsilon \tag{1.1}
$$

où E est le module d'Young, identifié comme la pente de la droite.

**Seconde phase** Elle est caractérisée par $\sigma > \sigma_0$. Le matériau gagne le domaine plastique. Pour observer l'irréversibilité liée au comportement plastique, il est nécessaire d'effectuer une décharge de l'éprouvette. La décharge effectuée à partir d'une contrainte $\sigma_1$ supérieure à la limite d'élasticité initiale $\sigma_0$ ne se fait pas suivant la courbe de charge mais suivant une droite de pente $E$ : c'est de nouveau un comportement élastique et réversible. On identifie donc une partition *additive* de la déformation uniaxiale en parties élastique et plastique :

$$
\epsilon = \epsilon^e + \epsilon^p \tag{1.2}
$$

où $\epsilon^e$ est la déformation élastique, réversible (récupérable à la décharge), liée à la contrainte
au travers de la loi d'élasticité ($\sigma = E\epsilon^e$) et où $\epsilon^p$ est la déformation plastique, irréversible
ou résiduelle.

On observe qu'il n'existe plus de relation *biunivoque* entre $\sigma$ et $\epsilon$ dû à l'irréversibilité de la déformation plastique et donc $\underline{\sigma} = f(\underline{\epsilon})$. Ce type de relation doit être remplacé par une loi de type « taux » traduisant les effets d'histoire associés à l'irréversibilité de la déformation plastique :

$$
\dot{\epsilon}^p = f(\sigma) \tag{1.3}
$$

Plus exactement, comme le temps intervenant ici n'est qu'un temps cinématique, on peut
réécrire cette même équation sous forme incrémentale $d\epsilon^p = f(\sigma)$ où $d\epsilon^p$ est un incrément
élémentaire de déformation plastique.

**Écrouissage** Si on recharge l'éprouvette après la décharge qui nous a permis d'identifier une déformation irréversible, la plasticité redémarre à $\sigma_1$ (effet lié à la réversibilité de la décharge-recharge) ; en conséquence, la limite d'élasticité est passée de $\sigma_0$ à $\sigma_1$ : on a *à écrouit* le matériau. L'écrouissage se traduit par le fait que la limite d'élasticité dépend de la déformation plastique.

Si la décharge se poursuit jusqu'à repasser en plasticité, cette fois-ci en compression, on
observe généralement qu'elle s'effectue à $\sigma_2 > -\sigma_0$. Aussi, l'écrouissage du matériau en
traction a modifié la limite d'élasticité en compression : c'est l'effet *Bauschinger* illustré sur
la figure 1.2. Cet effet n'intervient qu'après un premier chargement. La plasticité rompt la
symétrie traction-compression.
<figure><img src="image_13.png" /><figcaption>Figure 1.2 – Effet Bauschinger</figcaption></figure>

Comme évoqué dans l'introduction, ce type de comportement est principalement observé pour les matériaux métalliques, pour lesquels le tableau 1.1 liste quelques ordres de grandeur pour la limite d'élasticité en traction initiale et le module d'Young identifiés dans le cas d'alliages usuels en construction mécanique.

<table><caption>Tableau 1.1 – Propriétés de quelques matériaux métalliques</caption><thead><tr><td>Matériau</td><td>Module d'Young</td><td>Limité d'élasticité en traction</td></tr></thead><tbody><tr><td>Acier de construction non allié</td><td>210 000 MPa</td><td>285-355 MPa</td></tr><tr><td>Alliage d'aluminium</td><td>74 000 MPa</td><td>240 MPa</td></tr><tr><td>Acier faiblement allié à haute résistance</td><td>210 000 MPa</td><td>750-1450 MPa</td></tr><tr><td>Titane</td><td>116 000 MPa</td><td>260 MPa</td></tr><tr><td>Cuivre</td><td>100 000 MPa</td><td>40 MPa</td></tr></tbody></table>

### 1.1.2 Sollicitation multiaxiale

Bien que l'essai de traction permette d'avoir une première idée du matériau testé, il n'est malheureusement pas suffisant pour caractériser complètement le comportement de ce matériau. Les états de contraintes multidimensionnels auxquels peut être soumis un élément de matière dans une structure en service nécessite de maîtriser le comportement multiaxial du matériau, et ceci afin de préserver l'intégrité des structures conçues. Aussi, on utilise des essais qui permettent de générer ces sollicitations multiaxiales pour notamment identifier expérimentalement les surfaces seuil de plasticité dont nous parlerons en section 1.3.

La caractérisation multidimensionnelle du comportement relève rapidement du défi à mesure que l'état de contrainte se complique (dès que croît le nombre de composantes non nulles du tenseur de contrainte), aussi les essais les plus courants permettent d'identifier le comportement dans un plan donné; parmi les essais générant des états de contrainte multiaxiaux, on peut citer parmi les plus courants :

- l'essai de traction/compression-pression interne sur tube mince dans le plan $(\sigma_{\theta\theta}, \sigma_{zz})$, figure 1.3(a);

- l'essai de traction/compression-torsion d'un tube mince dans le plan $(\sigma, \tau)$, figure 1.3(b);

- l'essai de traction biaxiale sur éprouvette cruciforme dans le plan $(\sigma_1, \sigma_2)$, figure ??.

Ces essais tentent, dans la mesure de leur bon dimensionnement, de présenter au moins une zone utile dans laquelle l'état de contrainte est homogène pour faciliter le dépouillement de l'essai. Ces essais classiques ne suffisent pas toujours pour caractériser le comportement multiaxial d'un
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 1.3 –** Essais multi-axiaux</figcaption></figure>

matériau, et on est parfois amené à définir des essais ne présentant pas de zone dans laquelle l'état de contrainte est homogène, on parlera alors d'essais « structure ». Le dépouillement et l'interprétation de tels essais nécessitent la mise en œuvre d'analyses inverses et donc l'appel à un code de calcul, rendant leur coût particulièrement élevé.

Les sollicitations multi-axiales peuvent aussi être générées au travers de chargement proportionnels ou non-proportionnels. Pour ces derniers, les paramètres cinématiques de la machine d'essai sont pilotés de façon indépendante et non-proportionnelle. Ces chargements permettent, là encore, de mettre en évidence que l'état final de l'élément de matière dépend du trajet de chargement. Les effets d'histoire déterminent donc l'état final du matériau constitutif.

## 1.2 Notions sur les origines physiques

La plasticité observée à l'échelle macroscopique au travers de l'irréversibilité d'une partie de la déformation sur l'essai de traction est la conséquence de mécanismes de glissement à l'échelle microscopique. La plasticité est un phénomène intervenant dans les métaux, pour lesquels il existe une structure cristalline, c'est-à-dire que les atomes forment un réseau périodique dans un cristal. Suivant la construction du matériau, il peut être composé d'un seul grain (c'est alors un monocristal), par exemple pour les aubes de turbines, ou plus généralement fait d'un polycristal, chaque grain ayant son orientation propre et relié aux autres par des « joints de grain ».

Dans un cristal, il existe une position stable des atomes due à un minimum d'énergie potentielle. Autrement dit, il existe des forces inter-atomiques d'origine électromagnétique qui sont telles que ce minimum d'énergie est atteint au repos. En élasticité, la distorsion du réseau est réversible, ces forces inter-atomiques équilibrent le chargement imposé au matériau. Lorsque le chargement devient trop important, on assiste à un glissement des plans atomiques les uns sur les autres qui subsiste à la décharge : ce phénomène est illustré sur la figure 1.4. Cette vision du monocristal parfait explique :

<figure><img src="image_16.png" /><figcaption>**Figure 1.4 –** Charge et décharge d'un cristal</figcaption></figure>
— la réversibilité de l’élasticité : en effet celle-ci se traduit par une déformation et une distorsion du réseau, mais sans changer l’agencement relatif des atomes du cristal;

— l’irréversibilité de la plasticité : cette irréversibilité se traduit par un glissement des plans atomiques les uns sur les autres, et donc change « l’agencement relatif » des atomes dans le cristal;

— l’additivité des déformations ($\epsilon = \epsilon^e + \epsilon^P$) de par la structure du réseau.

Cependant, cette vision n’explique pas :

— l’écrouissage

— l’effet Bauschinger

Une autre difficulté subsiste : cette vision ne conduit pas à une estimation correcte de la limite d’élasticité du matériau. En effet, considérons, figure 1.5, deux plans atomiques soumis à un

<figure><img src="image_17.png" /><figcaption>**Figure 1.5** – Cisaillement d’un réseau cristallin</figcaption></figure>

cisaillement de sorte que le déplacement du plan supérieur atteigne la moitié de la distance inter-atomique dans la direction du cisaillement imposé, invitant à penser que cette position correspond à la limite avant de rompre les liaisons inter-atomiques entre les deux plans et d’en créer de nouvelles avec les atomes décalés d’un rang le long du plan cristallin.

Le glissement calculé $\gamma_{12} = 2\epsilon_{12}$ vaut $1/2$ (pas petit !). La contrainte de cisaillement associée s’écrit $\sigma_{12} = 2\mu\epsilon_{12} = \mu/2$, où $\mu$ est le second coefficient de Lamé, et fournit donc une estimation de l’ordre de grandeur de la limite d’élasticité $\sigma_0$. Pour un acier, $\mu$ est de l’ordre de 70 000 MPa ($E = 200\,000$ MPa et $v = 0,3$), donc $\sigma_0 \sim \mu/2 \sim 35\,000$ MPa; cette estimation est trop grande d’un facteur 100 ! Ainsi, l’hypothèse d’une séparation en bloc le long d’un plan atomique n’est pas réaliste.

L’explication tient en fait à l’existence de défauts au sein du cristal qu’on appelle des *dislocations* [1], illustrées sur la figure 1.6. Une dislocation est un point dans le plan mais une

<figure><img src="image_18.png" /><figcaption>**Figure 1.6** – Mouvement d’une dislocation le long d’un plan de glissement</figcaption></figure>

ligne dans un réseau tridimensionnel, et correspond donc à une ligne de lacunes dans le réseau atomique de base. Une dislocation va se propager de proche en proche jusqu’à déboucher, provoquant ainsi une « marche » analogue à celle que provoquerait un glissement en bloc d’un

<footnote>[1] SUQUET. *Rupture et plasticité*. École Polytechnique, 2004.</footnote>
<figure><img src="image_19.png" /></figure>

plan atomique. La contrainte nécessaire pour propager ce défaut est bien inférieure à celle nécessaire pour effectuer un mouvement en bloc. On peut aussi invoquer l'analogie du « pli dans le tapis »[1] illustrée sur la figure 1.7 : en créant un pli dans le tapis, puis en déplaçant le pli, les forces mises en jeu sont nettement moins importantes que celles nécessaires pour vaincre le frottement sur toute sa base, puisqu'il suffit seulement de vaincre le frottement sous une partie du tapis. L'effet final est le même lorsque le pli atteint l'extrémité du tapis.

<figure><img src="image_20.png" /><figcaption>Figure 1.7 – Analogie du « pli dans le tapis »</figcaption></figure>

Enfin, la présence de défauts explique le phénomène d'écrouissage : c'est le blocage des dislocations sur les obstacles. Plus on déforme le matériau, plus les dislocations se bloquent sur des obstacles. De plus, les dislocations reviennent plus facilement qu'elles ne continuent : ceci explique l'effet Baushinger.

D'autres conséquences ont lieu pour les métaux :

— La déformation plastique macroscopique est générée par des déformations de cisaillement à l'échelle locale, et donc sans variation de volume. Ceci entraîne la condition d'incompressibilité plastique :

$$ \mathrm{tr} \, \boldsymbol{\varepsilon}^{\mathrm{P}} = 0 \qquad (1.4) $$

— Les contraintes hydrostatiques (multiples de l'identité) n'ont aucun effet sur l'évolution des déformations plastiques du matériau.

## 1.3 Critères de plasticité

Revenons à présent à l'échelle de l'élément de matière, c'est-à-dire à l'échelle macroscopique. Pour définir quand les mécanismes de la plasticité sont activés, il faut introduire la notion de critère de plasticité. Dans le cas unidimensionnel de l'essai de traction, la surface de plasticité se réduit à deux points et le domaine d'élasticité est une droite :

— $|\sigma| < \sigma_0$ : élasticité

— $\sigma = \pm \sigma_0$ : plasticité possible

Dans le cas tridimensionnel, on peut se demander pour quelle « magnitude » du tenseur des contraintes $\sigma$ la plasticité apparaît. On introduit un critère portant sur une fonction de charge $f(\sigma)$ à valeur scalaire, définie de façon conventionnelle comme étant inférieure ou égale à zéro, de sorte que :

— $f(\sigma) < 0$ : élasticité

— $f(\sigma) = 0$ : plasticité possible

Cette fonction définit le domaine dans l'espace des contraintes à six dimensions (ou dans l'espace des contraintes principales à trois dimensions) à l'intérieur duquel toute variation de contrainte n'engendre que des variations de déformation élastique. L'ensemble des états de contrainte $\sigma$ tels que $f(\sigma) = 0$ définit la surface de charge. La fonction de charge permet de transformer un état de contrainte complexe à six dimensions en une valeur scalaire équivalente pouvant être comparée à une donnée expérimentale facile d'accès, la limite d'élasticité en traction.
<figure><img src="image_21.png" /></figure>

Cette fonction de charge est définie en fonction des caractéristiques propres du matériau considéré. Nous ne nous intéressons qu'aux matériaux métalliques, supposés abusivement vierges de toute déformation initiale, pour lesquels ces caractéristiques principales sont l'isotropie et une indépendance du comportement vis-à-vis de la contrainte hydrostatique :

— l'isotropie impose que la frontière du domaine d'élasticité soit invariante par changement de repère. Par conséquent, la fonction de charge $f(\sigma)$ ne dépend que des invariants du tenseur des contraintes $\sigma$ :

$$ 
\begin{align} 
I_1(\sigma) &= \operatorname{tr} \sigma \nonumber \\ 
I_2(\sigma) &= \frac{1}{2} [(\operatorname{tr} \sigma)^2 - \sigma : \sigma] \tag{1.5} \\ 
I_3(\sigma) &= \det \sigma \nonumber 
\end{align} 
$$

Le premier invariant qualifie, dans le critère, une sensibilité à la pression hydrostatique, et le second invariant, à l'énergie due au cisaillement. Enfin, si le troisième invariant n'intervient pas, la surface seuil est de révolution autour de l'axe hydrostatique;

— on observe généralement, dans les matériaux métalliques, l'incompressibilité plastique et une indépendance du comportement vis-à-vis de la contrainte hydrostatique. Cette indépendance s'écrit :

$$ f(\sigma_1 + \lambda, \sigma_2 + \lambda, \sigma_3 + \lambda) = f(\sigma_1, \sigma_2, \sigma_3), \quad \forall \lambda \qquad (1.6) $$

Par conséquent, la fonction $f$ est invariante par translation de vecteur directeur (1, 1, 1) et la surface de charge est un cylindre d'axe $\Delta(1, 1, 1)$. Il suffit alors de connaître une section de cette surface par un plan normal à l'axe directeur : ce plan est le plan du déviateur, on travaille alors avec les invariants du déviateur des contraintes $s$ qui s'écrit par définition :

$$ s = \sigma - \frac{1}{3}(\operatorname{tr} \sigma)\mathbf{1} \qquad (1.7) $$

Suivant le matériau considéré ou le procédé de mise en forme de la structure, l'isotropie et/ou l'indépendance du comportement vis-à-vis de la contrainte hydrostatique peuvent être remis en cause, en voici trois exemples :

— un matériau composite (stratifié par exemple) est par construction anisotrope et son comportement peut non seulement dépendre de la direction de chargement mais aussi du type de sollicitation traction/compression;

— le laminage d'une tôle, même à partir d'un matériau vierge pouvant être considéré comme isotrope, induit des directions préférentielles par déformation plastique et conduit à un comportement orthotrope qu'il faut considérer dans le critère de plasticité retenu pour le matériau constitutif de cette tôle;

— enfin, un matériau poreux présente une sensibilité à la pression hydrostatique dont il convient de rendre compte.

Les critères plus complexes anisotropes et/ou dépendant de la pression hydrostatique ne sont pas présentés dans ce document[2-4]. On présente, dans la suite, les critères usuels de von Mises et Tresca qui sont, de loin, les plus utilisés pour les matériaux métalliques.

<footnote>[2] SALENÇON. De l'élastoplasticité au calcul à la rupture. Éditions de l'École polytechnique, 2002.</footnote>

<footnote>[3] BESSON, CAILLETAUD, CHABOCHE et FOREST. Mécanique non linéaire des matériaux. Hermès, 2001.</footnote>

<footnote>[4] LEMAITRE, CHABOCHE, BENALLAL et DESMORAT. Mécanique des matériaux solides. Dunod, 2009.</footnote>
<figure><img src="image_22.png" /></figure>

### 1.3.1 Critère de von Mises

Le critère de von Mises considère que le seuil de plasticité est lié à l'énergie de déformation élastique de cisaillement. L'énergie de déformation élastique est définie à partir de la puissance des efforts intérieurs :

$$W^e = \int_{0}^{t} -P_{\text{int}} dt' = \int_{0}^{t} \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} dt' \quad \text{avec} \quad \dot{\boldsymbol{\epsilon}} = \dot{\boldsymbol{\epsilon}}^e \qquad (1.8)$$

En introduisant les parties déviatrices des tenseurs de contrainte et de déformation élastique dans l'expression de l'énergie de déformation élastique, on identifie la partie de cette énergie liée au cisaillement qu'on calcule pour un matériau élastique isotrope :

$$W_D^e = \int_0^t \mathbf{s} : \dot{\mathbf{s}} dt' = \frac{\mathbf{s} : \mathbf{s}}{4G} \qquad (1.9)$$

où $G$ est le module d'élasticité de cisaillement de Coulomb. Aussi le critère de Mises met en relation cette énergie de cisaillement avec une valeur seuil $W_y$ qui reste à identifier :

$$W_D^e \le W_y \qquad (1.10)$$

On établit généralement le lien entre cette énergie seuil et la limite d'élasticité en traction en identifiant ce critère sur l'essai de traction, l'état de contrainte et le déviateur s'écrivant pour cet essai de la façon suivante :

$$\boldsymbol{\sigma} = \sigma \mathbf{e}_1 \otimes \mathbf{e}_1 ; \quad \mathbf{s} = \frac{2}{3}\sigma \mathbf{e}_1 \otimes \mathbf{e}_1 - \frac{1}{3}\sigma \mathbf{e}_2 \otimes \mathbf{e}_2 - \frac{1}{3}\sigma \mathbf{e}_3 \otimes \mathbf{e}_3 \qquad (1.11)$$

le calcul donne :

$$\frac{\mathbf{s} : \mathbf{s}}{4G} = \frac{\sigma_y^2}{6G} \qquad (1.12)$$

La surface de charge $f(J_2, \sigma_y) = 0$ peut donc se réécrire en introduisant la contrainte équivalente $\sigma_{eq} = \sqrt{3J_2} = \sqrt{3s : s/2}$ (basée sur le second invariant du déviateur des contraintes $J_2 = s : s/2$) de la façon suivante :

$$f = \sigma_{\text{eq}} - \sigma_y = 0 \qquad (1.13)$$

Tout état de contrainte tridimensionnel vérifiant le critère (1.13) est équivalent au sens de von Mises à l'état unidimensionnel défini par le seuil de plasticité $\sigma_y$. Ce critère peut aussi s'exprimer dans l'espace des contraintes à six dimensions :

$$\frac{1}{2}[(\sigma_{11} - \sigma_{22})^2 + (\sigma_{22} - \sigma_{33})^2 + (\sigma_{33} - \sigma_{11})^2 + 6(\sigma_{12}^2 + \sigma_{23}^2 + \sigma_{13}^2)] - \sigma_y^2 = 0 \quad (1.14)$$

ou dans l'espace des contraintes principales à trois dimensions :

$$\frac{1}{\sqrt{2}}[(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_3 - \sigma_1)^2]^{1/2} = \sigma_y \qquad (1.15)$$

C'est l'équation d'un cylindre à base circulaire d'axe directeur $\Delta(1, 1, 1)$ et de rayon $\sqrt{2/3}\sigma_y$ illustré sur la figure 1.8.
<figure><img src="image_23.png" /><figcaption>Figure 1.8 – Critères de Mises et Tresca dans le plan des contraintes</figcaption></figure>

### 1.3.2 Critère de Tresca

L'idée du critère de Tresca, ou critère du cisaillement maximal, est que la plasticité n'est plus liée à l'énergie mais est provoquée par la contrainte de cisaillement maximale. Aussi le critère porte sur le maximum de la contrainte de cisaillement sur toutes les facettes possibles. À une facette donnée de normale **n**, le vecteur contrainte peut se décomposer de la façon suivante :

$$ T(\mathbf{n}) = \sigma \cdot \mathbf{n} = \sigma \mathbf{n} + \tau \qquad (1.16) $$

où $\sigma$ est la contrainte normale à la facette et $\tau$ le vecteur contrainte tangentielle évoluant dans le plan tangent. L'état de contrainte tridimensionnel en un point peut être représenté de façon plane au moyen des tricercles de Mohr schématisés sur la figure 1.9. Chaque cercle est le lieu géométrique des points décrits par l'extrémité du vecteur contrainte lorsque la normale **n** à la facette considérée varie dans un plan donné. Le tricercle de Mohr est construit à partir des contraintes principales $\sigma_1$, $\sigma_2$ et $\sigma_3$ données dans le repère principal ($e_1$, $e_2$, $e_3$). Le domaine coloré de la figure 1.9 est le lieu géométrique pouvant être décrit par l'extrémité du vecteur contrainte (ou état admissible des contraintes) lorsque la normale **n** à la facette considérée varie dans l'espace.

<figure><img src="image_24.png" /><figcaption>Figure 1.9 – Cercles de Mohr : $\mathbf{n} \in (e_1, e_2)$ [---], $\mathbf{n} \in (e_1, e_3)$ [---] et $\mathbf{n} \in (e_2, e_3)$ [---]</figcaption></figure>

Le critère de Tresca porte sur la contrainte de cisaillement maximale qui s'écrit comme la différence des contraintes principales extrêmes :

$$ \max \|\tau\| = \frac{1}{2} |\sigma_3 - \sigma_1| \qquad (1.17) $$
Le critère de Tresca s'écrit alors :

$$f(\boldsymbol{\sigma}) = \sup_{i,j} |\sigma_i - \sigma_j| - \sigma_0 \qquad (1.18)$$

De la même façon que le critère de von Mises, le critère de Tresca est invariant par ajout d'une contrainte hydrostatique puisque basé sur le cisaillement maximal. Comme indiqué sur la figure 1.8, il définit dans l'espace des contraintes principales un cylindre de vecteur directeur $\Delta(1, 1, 1)$ dont la projection dans le plan du déviateur est un hexagone inscrit dans le cylindre de von Mises.

### 1.3.3 Comparaison des critères de von Mises et Tresca

La figure 1.10 présente une comparaison des surfaces de charge liées aux critères de Mises et de Tresca dans le plan du déviateur. On note que les états admissibles en traction et compression pures dans les trois directions sont identiques pour les deux critères, seul les états mécaniques comprenant du cisaillement diffèrent. Un moyen simple pour comparer les deux critères consiste

<figure><img src="image_25.png" /><figcaption>**Figure 1.10** – Critères de von Mises [○] et Tresca [□] dans le plan du déviateur : traction pure [○], compression pure [●] et cisaillement pur [□]</figcaption></figure>

à considérer un état de cisaillement simple défini de la façon suivante :

$$\boldsymbol{\sigma} = \tau(\mathbf{e}_1 \otimes \mathbf{e}_2 + \mathbf{e}_2 \otimes \mathbf{e}_1) \qquad (1.19)$$

Calculons à présent pour quelle valeur de $\tau$ chacun des critères est atteint :

**Mises** Le calcul de la contrainte équivalente définie par (1.13) conduit à $\sqrt{3}\tau$. L'application du critère s'écrit $|\tau| = \sigma_y/\sqrt{3}$ et la limite d'élasticité en cisaillement simple calculée par Mises vaut $\sigma_y/\sqrt{3}$.

**Tresca** La contrainte équivalente de Tresca s'écrit comme la différence des contraintes principales extrémales $|\sigma_3 - \sigma_1|$, qui, pour un état de cisaillement simple, sont opposées et égales respectivement à $\tau$ et $-\tau$. L'application du critère conduit à $\tau = \sigma_y/2$ et la limite d'élasticité en cisaillement simple calculée par Tresca vaut $\sigma_y/2$ : elle est donc plus faible que celle de von Mises, c'est-à-dire plus restrictive<fnref n="1" /> du point de vue de la sécurité.

La figure 1.11 compare les deux critères dans les plans $(\sigma, \tau)$ et $(\sigma_1, \sigma_2)$. Ces représentations permettent d'identifier graphiquement les limites d'élasticité en cisaillement simple calculées par ces deux critères pour cet état de chargement.

<footnote n="1">¹ Conservative en anglais.</footnote>
<figure><img src="image_26.png" /></figure>

Les données expérimentales identifiées pour les surfaces de charge dépendent de l'alliage testé, certaines étant situées sur la frontière entre les deux critères. Les métaux très ductiles sont en général plus proches de von Mises. La figure 1.12 présente des expériences de traction-torsion menées sur du cuivre quasi-pur[5]. Il s'avère que le critère de von Mises prédit fort bien l'expérience.

<figure><img src="image_27.png" /><figcaption>**Figure 1.12** – Essai de traction-torsion sur cuivre : surface seuil initiale. En trait plein la surface seuil prédite par le critère de von Mises</figcaption></figure>

Dans les applications, le critère de von Mises est généralement préféré à celui de Tresca en raison de la régularité de sa surface de charge qui permet une plus grande facilité d'emploi dans les calculs numériques (traités au quatrième chapitre). Le critère de Tresca est plus délicat d'utilisation en raison des points singuliers de sa surface de charge qui se comporte comme des points attracteurs dans les algorithmes numériques classiques d'intégration de lois élastoplastiques.

## 1.4 Modélisation de l'écrouissage

Lorsque l'état de contrainte atteint la surface de charge initiale et que le chargement du matériau se poursuit, celui-ci se déforme plastiquement. Selon le matériau, cet état de charge peut progresser dans l'espace des contraintes en modifiant la géométrie de la surface de charge initiale, de sorte que cet état de charge reste toujours sur la surface de charge. On parle d'écrouissage au sens où la charge a modifié la géométrie du domaine élastique. Une décharge

<footnote>[5] BUI. « Évolution de la frontière du domaine élastique des métaux avec l'écrouissage plastique et le comportement élastoplastique d'un agrégat de cristaux cubiques ». Mémoire de l'Artillerie Française 1 (1971), pages 141-165.</footnote>
de l'éprouvette montrera d'ailleurs un retour dans ce domaine élastique qui vient d'être modifié. L'écrouissage dépend bien entendu du trajet de charge parcouru hors du domaine d'élasticité initial. Les modifications subies par la surface de charge peuvent être de trois sortes :

- une expansion isotrope de la surface;

- une translation du centre du domaine d'élasticité initial. Ce déplacement est à l'origine de l'effet Bauschinger observé sur un essai unidimensionnel : cette translation rompt la symétrie traction/compression du seuil de plasticité;

- une distorsion ou un changement de forme de la surface. Initialement elliptique, comme illustré sur la figure 1.12 dans un plan $(\sigma, \tau)$, cette surface peut se distordre jusqu'à présenter un point plus ou moins anguleux dans la direction de chargement.

Pour d'autres matériaux, la nouvelle surface seuil obtenue après écoulement plastique s'avère être identique à la surface de charge initiale, on parlera de matériau élastique *parfaitement plastique*.

Dans la suite, nous nous limitons à la modélisation du changement de taille de la surface (écrouissage isotrope) et à la translation du centre du domaine (écrouissage cinématique).

### 1.4.1 Écrouissage isotrope

Ce type d'écrouissage tient compte d'une expansion du domaine d'élasticité avec la déformation plastique. Cette expansion correspond à une homothétie de centre $0$. Aussi, on ne tient pas compte de l'effet Bauschinger, la surface reste centrée sur l'origine. Supposons que la surface de charge puisse s'écrire sous la forme :

$$f(\sigma, p) = \sigma_{\text{eq}} - \sigma(p) \qquad (1.20)$$

où $\sigma(p)$ est une fonction croissante de la déformation plastique équivalente cumulée, cette dernière étant définie de la façon suivante :

$$p = \int_{0}^{t} \dot{p}(\tau) d\tau; \quad \dot{p} = \sqrt{\frac{2}{3} \dot{\varepsilon}^{\text{P}} : \dot{\varepsilon}^{\text{P}}} \qquad (1.21)$$

La fonction $\sigma(p)$ est supposée connue expérimentalement, on verra qu'on peut en proposer plusieurs modèles. Cette fonction scalaire est astreinte à valoir la limite d'élasticité initiale lorsque la déformation plastique équivalente cumulée est nulle, soit formellement :

$$\sigma(p = 0) = \sigma_0 \qquad (1.22)$$

La figure 1.13 montre l'expansion du domaine d'élasticité décrite par ce type d'écrouissage.

<figure><img src="image_28.png" /><figcaption>Figure 1.13 – Expansion du domaine d'élasticité de von Mises : écrouissage isotrope</figcaption></figure>
L'écriture la plus simple pouvant être considérée $\sigma(p)$ est une fonction affine de la déformation plastique cumulée $p$, de sorte que $\sigma(p) = Qp + \sigma_0$, Q étant une constante. On peut aussi écrire cette fonction comme une loi puissance $\sigma(p) = Kp^n + \sigma_0$, c'est la loi de Ramberg-Osgood.

## 1.4.2 Écrouissage cinématique

L'écrouissage cinématique correspond à une translation du centre du domaine d'élasticité, sans changement de forme ni de taille, et permet donc de prendre en compte l'effet Bauschinger identifié lors de l'essai unidimensionnel. On introduit alors le tenseur déviateur $\mathbf{X}$ (tr $\mathbf{X} = 0$) pour décrire la position du centre du domaine d'élasticité. Ce tenseur est nul pour le domaine initial centré en $\mathbf{0}$. Aussi le critère actuel de plasticité s'applique à présent à l'état de contrainte « décentré » $\sigma - \mathbf{X}$, et s'écrit :

$$f(\sigma - \mathbf{X}) \le 0 \qquad (1.23)$$

Par exemple la contrainte équivalente de von Mises prend la forme suivante :

$$\sigma_{\text{eq}}^{\text{VM}} = \sqrt{\frac{3}{2}(\mathbf{s} - \mathbf{X}) : (\mathbf{s} - \mathbf{X})} \qquad (1.24)$$

où on rappelle que $\mathbf{s}$ est le déviateur des contraintes. Physiquement, $\mathbf{X}$ correspond à l'existence de contraintes internes à l'échelle microscopique dues aux défauts et aux hétérogénéités qui représentent des obstacles à la migration des dislocations dans les grains.

<figure><img src="image_29.png" /><figcaption>**Figure 1.14** – Translation du domaine d'élasticité de von Mises : écrouissage cinématique</figcaption></figure>

La figure 1.14 illustre la translation du domaine d'élasticité paramétrée par le tenseur $\mathbf{X}$. Il existe de nombreux modèles d'évolution de $\mathbf{X}$ dont certains assez complexes dont nous reparlerons au second chapitre. Le modèle le plus simple consiste à piloter son évolution de façon linéaire avec la vitesse de déformation plastique :

$$\dot{\mathbf{X}} = H \dot{\mathbf{e}}^p \qquad (1.25)$$

où $H$ est un paramètre. L'évolution linéaire est connue sous le nom de loi de Prager (1950).

**Essai de traction** Dans le cas de l'essai de traction, le tenseur des contraintes $\sigma$ et $\mathbf{X}$ s'écrivent respectivement :

$$\sigma = \sigma \mathbf{e}_1 \otimes \mathbf{e}_1 \quad ; \quad \mathbf{X} = \frac{X}{3}(2\mathbf{e}_1 \otimes \mathbf{e}_1 - \mathbf{e}_2 \otimes \mathbf{e}_2 - \mathbf{e}_3 \otimes \mathbf{e}_3) \qquad (1.26)$$

Le calcul de la contrainte équivalente de von Mises (1.24) donne :

$$\sigma_{\text{eq}}^{\text{VM}} = |\sigma - \mathbf{X}| \qquad (1.27)$$
<figure><img src="image_30.png" /><figcaption>**Figure 1.15 – Écrouissage cinématique : essai de traction**</figcaption></figure>

soit, par application du critère, l'état de contrainte actuel s'exprime comme $\sigma = \sigma_0 + X = \sigma_0 + H \varepsilon^P$. Ce décentrement correspondant à l'effet Bauschinger est illustré dans le cas de l'essai de traction sur la figure 1.15.

L'expérience montre que ce type d'écrouissage est de première importance dès qu'une structure subit des chargements cycliques, dans lesquels des états de contrainte de traction/compression peuvent se succéder dans chaque cycle. Si ces chargements sont suffisamment importants pour générer un écoulement plastique (par exemple en fatigue oligocyclique), la prise en compte d'une dissymétrisation traction/compression de la surface d'élasticité au travers de sa translation s'avère nécessaire pour pouvoir prédire correctement le comportement de l'élément de matière. Dans la pratique, on utilise des modèles qui combinent les écrouissages isotrope et cinématiques pour de nombreux matériaux métalliques.

## 1.5 Règle d'écoulement plastique

La définition du comportement plastique se compose schématiquement de la réponse à deux questions :

1. Quand y-a-t-il déformation plastique ? La réponse est donnée par les critères faisant intervenir les fonctions de charge, objet de la section précédente.

2. Comment s'effectue la déformation plastique ? C'est la **règle d'écoulement plastique** que nous allons aborder dans cette section.

### 1.5.1 Partition de la déformation

S'appuyant sur les observations réalisées sur l'expérience de traction, figure 1.1, on *postule* dans le cas des petites déformations une partition additive de la déformation totale en parties élastique et plastique pour le cas tridimensionnel de façon analogue à l'observation unidimensionnelle :

$$ \varepsilon = \varepsilon^e + \varepsilon^P \qquad (1.28) $$

La déformation élastique $\varepsilon^e$ est réversible et donc récupérable à la décharge, la déformation plastique $\varepsilon^P$ est irréversible et non récupérable à la décharge.

**Remarque** Il n'existe pas de partition du déplacement analogue à celle postulée pour la déformation :

$$ \mathbf{u} = \mathbf{u}^e + \mathbf{u}^P : \text{FAUX} \qquad (1.29) $$

De façon plus précise, seule la déformation totale vérifie les équations de compatibilité, et donc est relié au déplacement par $\varepsilon = 1/2(\nabla \mathbf{u} + \nabla^T \mathbf{u})$. Les déformations élastique et plastique ne
vérifient pas les équations de compatibilité, on verra par la suite que ceci implique l'existence de contraintes résiduelles. Il en est de même pour les déformations thermiques en thermoélasticité, pouvant précisément conduire à l'existence de contraintes résiduelles d'origine thermique.

**Déformation élastique** La déformation élastique est reliée au tenseur des contraintes au travers de la loi de Hooke :

$$ \boldsymbol{\epsilon}^e = \boldsymbol{S} : \boldsymbol{\sigma} \qquad (1.30) $$

où $\boldsymbol{S}$ est le tenseur des souplesses (tenseur du quatrième ordre).

**Déformation plastique** Il nous faut à présent préciser la loi de comportement de la déformation plastique. On a cependant observé à la section 1.1.1 sur l'essai unidimensionnel qu'il n'existe plus de relation biunivoque entre la contrainte et la déformation dès qu'il y a écoulement plastique, il n'existe pas unicité de la déformation plastique à un même niveau de contrainte. Aussi, nous ne sommes désormais plus en mesure que de préciser un sens de l'évolution des déformations plastiques à défaut d'une loi biunivoque. En revanche, si on connaît non seulement la contrainte actuelle $\boldsymbol{\sigma}$ mais aussi toute l'histoire du chargement à partir de l'état initial, alors la donnée d'une variation de contrainte $d\boldsymbol{\sigma}$ suffit à déterminer la variation correspondante de la déformation $d\boldsymbol{\epsilon}$. C'est donc une loi de comportement incrémentale qui régit le comportement élastoplastique.

**Conditions de charge-décharge** L'écoulement plastique ne peut se produire que si le point représentatif de l'état de contrainte $\boldsymbol{\sigma}$ atteint la surface de charge (le critère est atteint) et y reste, de sorte que ce point ne puisse sortir de la surface ($f > 0$ est impossible). Cette dernière condition est appelée condition de cohérence (ou de consistance) et s'écrit :

$$ \dot{f} = 0 \quad \text{ou} \quad df = 0 \qquad (1.31) $$

ce qui implique que le point représentatif de l'état $\boldsymbol{\sigma} + d\boldsymbol{\sigma}$ reste sur la surface de charge. En résumé :

$$ f < 0 \quad \rightarrow \quad \text{comportement élastique} $$

$$ f = 0 \text{ et } \dot{f} = 0 \quad \rightarrow \quad \text{écoulement plastique } (\dot{\boldsymbol{\epsilon}}^P \neq 0) $$

$$ f = 0 \text{ et } \dot{f} < 0 \quad \rightarrow \quad \text{décharge élastique} $$

### 1.5.2 Principe de Hill

Le principe de Hill[6], ou principe du « travail » plastique maximal énonce la règle d'évolution de la vitesse de déformation plastique : considérons un état de contrainte $\boldsymbol{\sigma}$ provoquant un taux de déformation plastique $\dot{\boldsymbol{\epsilon}}^P$ et un état de contrainte $\boldsymbol{\sigma}^*$ plastiquement admissible, c'est-à-dire respectant le critère $f(\boldsymbol{\sigma}^*) \le 0$, alors :

$$ (\boldsymbol{\sigma} - \boldsymbol{\sigma}^*) : \dot{\boldsymbol{\epsilon}}^P \ge 0 \qquad (1.32) $$

En réécrivant cette expression $\boldsymbol{\sigma}^* : \dot{\boldsymbol{\epsilon}}^P \le \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^P$, on observe qu'à $\dot{\boldsymbol{\epsilon}}^P$ fixé, la « puissance plastique » $\boldsymbol{\sigma}^* : \dot{\boldsymbol{\epsilon}}^P$ est maximale pour $\boldsymbol{\sigma}^* = \boldsymbol{\sigma}$, qui est l'état de contrainte provoquant $\dot{\boldsymbol{\epsilon}}^P$. Autrement dit, le principe de Hill résulte de la maximisation de la puissance plastique, de sorte que :

$$ \mathcal{D}(\boldsymbol{\sigma}) = \max_{\boldsymbol{\sigma}^* \in \mathbb{P}} (\boldsymbol{\sigma}^* : \dot{\boldsymbol{\epsilon}}^P), \quad \mathbb{P} = \{\boldsymbol{\sigma} \in \mathcal{S}^{\text{ad}}, f(\boldsymbol{\sigma}) = 0\} \qquad (1.33) $$

<footnote>[6] HILL. *The Mathematical Theory of Plasticity*. Clarendon Press, 1967.</footnote>
Ce principe s'applique au cas des matériaux parfaitement plastiques. On observe d'autre part qu'en prenant $\boldsymbol{\sigma}^* = \mathbf{0}$, on tombe sur l'expression du second principe de la thermodynamique [section 2.1.2], aussi ce dernier s'avère moins restrictif que le principe de Hill.

**Conséquences du principe de Hill** Le taux de déformation plastique $\dot{\boldsymbol{\epsilon}}^P$ est nécessairement normal à la surface de charge, comme indiqué sur la figure 1.16 : c'est la loi de normalité. La

<figure><img src="image_31.png" /><figcaption>Figure 1.16 – Normalité de l'écoulement plastique</figcaption></figure>

normalité se traduit par l'expression suivante :

$$ \dot{\boldsymbol{\epsilon}}^P = \dot{\boldsymbol{\lambda}} \frac{\partial f}{\partial \boldsymbol{\sigma}} ; \quad \dot{\boldsymbol{\lambda}} \ge 0 \qquad (1.34) $$

où $\dot{\boldsymbol{\lambda}}$, dénommé le multiplicateur plastique, est positif ou nul. On observe aussi que la surface de charge $f$ joue le rôle d'un pseudo-potentiel plastique. La condition de normalité peut se déduire directement de la maximisation (1.33), mais on peut aussi l'appréhender en appliquant le principe de Hill (1.32) à partir d'un point $\boldsymbol{\sigma}$ de la surface de charge, et en choisissant un point $\boldsymbol{\sigma}^*$ infiniment proche, aussi sur la surface de charge. On peut alors déduire $\boldsymbol{\sigma}^*$ de $\boldsymbol{\sigma}$ par un tenseur $\mathbf{t}^*$ appartenant au plan tangent de la surface en $\boldsymbol{\sigma}$, de sorte que $\boldsymbol{\sigma}^* = \boldsymbol{\sigma} + k\mathbf{t}^*$, avec $k > 0$. Le même raisonnement peut être conduit avec $\boldsymbol{\sigma}^* = \boldsymbol{\sigma} - k\mathbf{t}^*$ comme point de départ, conduisant aux inégalités suivantes :

$$ k\mathbf{t}^* : \dot{\boldsymbol{\epsilon}}^P \ge 0 \quad \text{et} \quad -k\mathbf{t}^* : \dot{\boldsymbol{\epsilon}}^P \ge 0 \qquad (1.35) $$

si bien que

$$ \mathbf{t}^* : \dot{\boldsymbol{\epsilon}}^P = 0 \qquad (1.36) $$

cette dernière expression montrant bien que la vitesse de déformation plastique est portée par la normale à la surface de charge. Tout se ramène donc à présent à calculer le multiplicateur plastique $\dot{\boldsymbol{\lambda}}$, définissant l'intensité de l'écoulement.

Le domaine d'élasticité est *convexe*. Si la surface de charge n'est pas *convexe*, le domaine n'est pas entièrement du même côté de l'hyperplan tangent, le principe de Hill est violé. La figure 1.17 illustre la violation du principe (1.32) pour une surface non *convexe*. Ces conséquences impliquent un ensemble de conditions sur $\dot{\boldsymbol{\lambda}}$ et $f$, en effet :

— si $f < 0$, il y a élasticité, donc $\dot{\boldsymbol{\lambda}} = 0$ (pas d'écoulement plastique).

— si $f = 0$, il peut y avoir écoulement plastique, donc $\dot{\boldsymbol{\lambda}} \ge 0$. De façon plus précise, et en introduisant la condition de charge plastique :

— si $\dot{f} < 0$, il y a décharge élastique, donc $\dot{\boldsymbol{\lambda}} = 0$.

— si $\dot{f} = 0$, il y a charge plastique $\dot{\boldsymbol{\lambda}} \ge 0$.
<figure><img src="image_32.png" /><figcaption>**Figure 1.17** – Non convexité de la surface de charge</figcaption></figure>

Ces conditions peuvent se résumer par les conditions de complémentarité de Kühn-Tucker :

$$ \dot{\lambda} \ge 0 \quad ; \quad f \le 0 \quad ; \quad \dot{\lambda} f = 0 \qquad (1.37) $$

En résumé, le principe du travail plastique maximal est une propriété d'extrémalité caractéristique de la dissipation plastique. Il implique, outre la convexité du domaine actuel d'élasticité, la normalité du taux de déformation plastique vis-à-vis de ce domaine. Le principe du travail plastique maximal est une hypothèse qui est ou n'est pas vérifiée selon les matériaux. Elle est vérifiée en première approximation pour les métaux, elle n'est pas vérifiée pour les sols. Ce principe permet d'engendrer une classe de modèles : les matériaux standards ; ceux-ci permettent de traiter de nombreux problèmes (nous en reparlerons au second chapitre).

### 1.5.3 Calcul du multiplicateur plastique

Le multiplicateur plastique $\dot{\lambda}$ introduit dans la loi de normalité (1.34) est laissé indéterminé. Dans le cas d'un matériau à écrouissage positif, cette indétermination peut être levée par la condition de cohérence (1.31). Pour un matériau élastique parfaitement plastique, il est déterminé par la résolution du problème aux limites. En ce sens, on peut en faire l'analogie dans ce dernier cas avec la pression hydrostatique dans le cas d'un écoulement visqueux incompressible, celle-ci étant aussi indéterminée dans la loi de comportement mais pouvant être déterminée par la résolution du problème aux limites.

**Cas de l'écrouissage isotrope** On considère pour la suite le critère de von Mises (1.20). La loi de normalité (1.34) appliquée avec ce critère s'écrit :

$$ \dot{\varepsilon}^P = \frac{3}{2} \dot{\lambda} \frac{s}{\sigma_{eq}} \qquad (1.38) $$

Comme évoqué ci-dessus, la condition de cohérence permet dans le cas d'un matériau à écrouissage positif de lever l'indétermination du multiplicateur plastique; cette condition se développe de la façon suivante :

$$ 0 = \dot{f}(\sigma, p) = \frac{\partial f}{\partial \sigma} : \dot{\sigma} + \frac{\partial f}{\partial p} \dot{p} \qquad (1.39) $$

D'autre part, en combinant (1.38) et (1.21), il est possible de montrer que :

$$ \dot{p} = \dot{\lambda} \qquad (1.40) $$

En injectant (1.40) dans la condition de cohérence (1.39), il vient l'expression du multiplicateur plastique pour un tel matériau :

$$ \dot{\lambda} = \frac{\dot{\sigma}_{eq}}{\frac{d\sigma}{dp}(p)} \qquad (1.41) $$

La loi d'écoulement plastique pour le critère de von Mises avec écrouissage isotrope se résume donc à :

$$ \dot{\varepsilon}^P = \frac{3}{2} \dot{\lambda} \frac{s}{\sigma_{eq}} \quad ; \quad \dot{\lambda} = \frac{\dot{\sigma}_{eq}}{\frac{d\sigma}{dp}(p)} \qquad (1.42) $$
**Cas de l'écrouissage cinématique** L'expression de la contrainte équivalente en présence d'écrouissage cinématique (1.24) conduit à l'expression suivante de la loi de normalité :

$$ \dot{\boldsymbol{e}}^p = \frac{3}{2} \dot{\boldsymbol{\lambda}} \frac{\boldsymbol{s} - \boldsymbol{X}}{\sigma_{eq}(\boldsymbol{s} - \boldsymbol{X})} \qquad (1.43) $$

De la même manière, l'équation de cohérence permet de calculer le multiplicateur plastique, celle-ci se développe :

$$ 0 = \dot{f}(\boldsymbol{\sigma}, \boldsymbol{X}) = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + \frac{\partial f}{\partial \boldsymbol{X}} : \dot{\boldsymbol{X}} \qquad (1.44) $$

Soit en injectant la loi d'évolution de l'écrouissage (1.25), et la loi d'écoulement (1.43) dans (1.44), il vient l'expression du multiplicateur, la loi d'écoulement plastique pour le critère de von Mises avec écrouissage cinématique se résume donc à :

$$ \dot{\boldsymbol{e}}^p = \frac{3}{2} \dot{\boldsymbol{\lambda}} \frac{\boldsymbol{s} - \boldsymbol{X}}{\sigma_{eq}(\boldsymbol{s} - \boldsymbol{X})} ; \quad \dot{\boldsymbol{\lambda}} = \frac{(\boldsymbol{s} - \boldsymbol{X}) : \dot{\boldsymbol{\sigma}}}{H\sigma_{eq}} \qquad (1.45) $$

**Cas de la plasticité parfaite** Le cas de la plasticité parfaite peut être considéré comme un cas limite des cas précédents pour une pente d'écrouissage tendant vers zéro. Partant du cas de l'écrouissage isotrope (1.42), on a en plasticité parfaite :

$$ \frac{d\boldsymbol{\sigma}}{dp} = 0 \quad \text{mais aussi} \quad \dot{\boldsymbol{\sigma}}_{eq} = 0 \qquad (1.46) $$

il s'ensuit que le multiplicateur plastique $\dot{\boldsymbol{\lambda}}$ est indéterminé. Partant du cas de l'écrouissage cinématique (1.45), il vient en plasticité parfaite $H = 0$ et $\dot{\boldsymbol{\sigma}} = \boldsymbol{0}$ : on parvient donc à la même conclusion. Aussi on peut poser en plasticité parfaite que la vitesse d'écoulement plastique est colinéaire au déviateur des contraintes :

$$ \dot{\boldsymbol{e}}^p = \dot{\boldsymbol{\lambda}}\boldsymbol{s} ; \quad \dot{\boldsymbol{\lambda}} \ge 0 \qquad (1.47) $$

Le multiplicateur $\dot{\boldsymbol{\lambda}}$ sera déterminé par la résolution du problème aux limites.

### 1.5.4 Postulat de Drucker-Ilyushin et principe de Hill

Il existe une grande latitude possible pour la construction des lois d'évolution. La part d'arbitraire est généralement réduite par certains principes généraux parmi lesquels le principe d'objectivité (indépendance à l'observateur) et le deuxième principe de la thermodynamique. Cependant ces principes laissent encore un choix très large. Pour le réduire, un autre postulat fondé sur une inégalité et dit « quasi-thermodynamique »[7] peut être appliqué à l'élastoplasticité. Ce postulat dit que « le travail de déformation doit être non négatif dans tout cycle en déformation compatible avec la loi d'évolution » et s'écrit :

$$ W = \oint \boldsymbol{\sigma} : \dot{\boldsymbol{e}} dt \ge 0 \qquad (1.48) $$

Cela revient à dire qu'il faut fournir de l'énergie au matériau pour lui faire faire un cycle en déformation. Dans le cas d'un matériau élastique parfaitement plastique (pour lequel le principe de Hill s'applique), il est possible de montrer l'équivalence entre le postulat de Drucker-Ilyushin et le principe de Hill[7].

<footnote>[7] MARIGO. Plasticité et rupture. [cel-01374813]. École Polytechnique, 2011.</footnote>
<figure><img src="image_33.png" /><figcaption>**Figure 1.18** – Écrouissage positif ou nul impliqué par le postulat de Drucker-Ilyushin</figcaption></figure>

Le postulat de Drucker-Ilyushin a cependant une portée plus générale, et peut s'appliquer aux matériaux écrouissables. Il permet alors de restreindre l'éventail de choix des lois d'évolutions admissibles pour un matériau donné. Une autre façon d'écrire le postulat de Drucker[2] vient en recherchant le sup($\dot{\sigma} : \dot{\epsilon}^P$), $\forall \dot{\sigma}$. Dans le cas général, on peut écrire :

$$ \dot{\sigma} : \dot{\epsilon}^P \ge 0 \qquad (1.49) $$

Pour un matériau élastique parfaitement plastique, l'inégalité (1.49) devient une égalité. Une manière de l'interpréter est de considérer une courbe de traction, illustrée sur la figure 1.18, on observe alors que ce postulat interdit le phénomène d'écrouissage négatif ou adoucissement.
# Thermodynamique des milieux continus

L'expérience montre que les transformations internes des corps sont en partie irréversibles et s'accompagnent d'échanges thermiques et de variations de température. Il en résulte qu'une description satisfaisante des lois de comportement doit s'inscrire dans le cadre de la thermodynamique macroscopique des processus irréversibles. La déformation d'un milieu continu est une transformation thermodynamique particulière, elle doit donc s'opérer en accord avec les deux principes de la thermodynamique, qui fixe des limites aux échanges de chaleur et d'énergie mécanique des systèmes physiques. La thermodynamique des milieux continus à variables internes apparaît donc comme une démarche rigoureuse pour construire un modèle de comportement.

La plasticité représente à cet égard *un type de comportement s'accompagnant de phénomènes irréversibles*. L'objet de ce chapitre est d'intégrer au formalisme thermodynamique le modèle du solide élastoplastique présenté au premier chapitre sur la base de principes. Dans la section 2.1, on rappelle brièvement les principes de la thermodynamique. Dans la section 2.2, on s'appuie alors sur ces bases pour dégager la formulation des lois élastoplastiques en *lois d'état* et *lois complémentaires*. On montrera que pour une classe de matériaux appellés *matériaux standards généralisés*, chacune de ces lois dérive d'un potentiel, qu'on particularise au cas de la plasticité associée dans la section 2.3. Enfin, dans la section 2.4, on montrera que pour tenir compte de phénomènes plus complexes tels que l'écrouissage cinématique non linéaire, on est amené à devoir généraliser l'approche précédente à la plasticité *non associée*.

## 2.1 Rappels de thermodynamique

On s'intéresse ici à des systèmes fermés, en équilibre mécanique et thermique, et décrits par un certain nombre de variables d'état. On suppose aussi que les évolutions lors des transformations sont suffisamment lentes pour que les notions introduites puissent s'étendre à des systèmes hors équilibre.

### 2.1.1 Premier principe de la thermodynamique

L'énergie est *quelque chose* qui jauge la capacité d'un système à effectuer des transformations comme fournir du travail, modifier un mouvement, modifier la température, déformer un corps, changer l'état de la matière, etc. Elle possède une propriété fondamentale : elle se conserve. C'est ce que traduit le premier principe de la thermodynamique :

$$ \frac{d\mathcal{E}}{dt} = \mathcal{P}_{\text{ext}} + \dot{\mathcal{Q}} \qquad (2.1) $$

où $\mathcal{E}$ dénote l'énergie totale du système, c'est-à-dire l'énergie interne $E$ plus l'énergie cinétique $E_c$, $\mathcal{P}_{\text{ext}}$ est la puissance des efforts extérieurs exercés sur le système et $\dot{\mathcal{Q}}$ exprime le taux de
chaleur fournie au système.

L'énergie interne $E$ est une quantité introduite par ce premier principe. C'est une grandeur extensive, définie à une constante additive près, dont on peut connaître la variation entre deux instants. Elle correspond à de l'énergie stockée à l'intérieur du système<fnref n="1" />. On peut avoir une idée de la nature de cette énergie en considérant une transformation adiabatique, c'est-à-dire en fournissant un travail au système : ce travail tend alors à augmenter cette énergie par augmentation de l'agitation du milieu à l'échelle microscopique.

Pour un milieu continu, chacune de ces grandeurs est définie comme la sommation de quantités locales :

$$E = \int_{\Omega} \rho e \, d\Omega \qquad E_c = \int_{\Omega} \frac{1}{2} \rho \mathbf{v}^2 \, d\Omega \qquad (2.2)$$

$$P_{\text{ext}} = \int_{\partial\Omega} (\boldsymbol{\sigma} \cdot \mathbf{n}) \cdot \mathbf{v} \, dS + \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\Omega \qquad \dot{Q} = \int_{\partial\Omega} -\mathbf{q} \cdot \mathbf{n} \, dS + \int_{\Omega} r \, d\Omega \qquad (2.3)$$

La puissance des efforts extérieurs $P_{\text{ext}}$ est générée par des efforts appliqués sur le bord du domaine et par des forces volumiques appliquées au système. Le taux de chaleur $\dot{Q}$ est composé d'un premier terme associé à la conduction opérée par le bord et d'un second associé à des actions d'origine extérieure<fnref n="2" />.

Cependant, le premier principe (2.1) peut aussi s'écrire comme un bilan d'énergie interne, en introduisant la puissance des efforts intérieurs dans l'équation. Pour cela, on combine l'équation (2.1) au principe des puissances virtuelles appliqué au mouvement réel, c'est-à-dire le théorème de l'énergie cinétique, il vient alors :

$$\frac{dE}{dt} = -P_{\text{int}} + \dot{Q} \qquad (2.4)$$

où la puissance des efforts intérieurs $P_{\text{int}}$ est telle que :

$$P_{\text{int}} = - \int_{\Omega} \boldsymbol{\sigma} : \mathbf{D} \, d\Omega \qquad (2.5)$$

et $\mathbf{D}$ étant le taux de déformation eulérien. L'équation (2.4) s'écrit sous forme locale de la façon suivante :

$$\rho \dot{e} = \boldsymbol{\sigma} : \mathbf{D} - \operatorname{div} \mathbf{q} + r \qquad (2.6)$$

Cette équation de conservation de l'énergie exprime l'origine des variations de l'énergie interne du système par un apport au travers d'une puissance mécanique des efforts intérieurs et d'une puissance calorifique fournies au système.

## 2.1.2 Second principe de la thermodynamique

Le second principe de la thermodynamique établit la distinction entre transformations réversibles et irréversibles. Si le premier principe est une loi de bilan, comme le bilan de la quantité de mouvement par exemple, le second principe exprime surtout un principe d'évolution qui permet de prévoir le sens de l'évolution d'un système. Le second principe introduit aussi les notions de température et d'entropie. La température est liée à l'agitation des atomes à l'échelle microscopique, c'est une grandeur intensive, mesurée en Kelvin, et n'est jamais négative. L'entropie $S$ est une fonction d'état extensive qui ne peut que croître au cours d'une transformation réelle. D'une façon générale, l'entropie est une grandeur qui caractérise la capacité d'un système

<footnote n="1">¹ micro-mécanismes, agitation thermique, interaction atomique, etc.</footnote>

<footnote n="2">² effet Joule, réactions chimiques, etc.</footnote>
<figure><img src="image_34.png" /></figure>

physique à subir des transformations spontanées : plus grande est la valeur de l'entropie, plus faible est la capacité du système à se transformer. En fait, un système perd de sa capacité à évoluer à mesure qu'il évolue (son entropie croît). L'entropie est liée à ce que certains appellent la « qualité » de l'énergie disponible au sein d'un système. Cette « qualité » décroît à mesure que le système se transforme, dans le sens où l'énergie disponible devient de moins en moins utilisable facilement<fnref n="3" />. La mécanique statistique l'assimile quant à elle à la notion de désordre à l'échelle microscopique.

Considérons par exemple une transformation isotherme s'effectuant à la température T₀, passant d'un état E₁ à un état E₂, celle-ci vérifie alors :

$$S_2 - S_1 \geq \frac{Q}{T_0} \qquad (2.7)$$

où Q exprime la chaleur reçue par le système au cours du processus. Il est intéressant d'observer que, pour un *processus réversible*, l'inégalité (2.7) devient une égalité. L'énergie calorifique peut s'écrire $\delta Q = T_0 dS$, c'est-à-dire comme le produit d'une variable intensive — la température — et de la différentielle d'une variable extensive — l'entropie. Cette écriture, dans le cas d'un *processus réversible*, peut être rapprochée à d'autres formes d'énergie, comme par exemple le travail d'un ressort ($\delta W = F d\ell$). Aussi, pour ce type de processus et dans le cas de cette écriture, l'entropie apparaît comme une variable « flux » activée par une « force thermodynamique », ici la température.

Dans le cas général des transformations *irréversibles*, il y a création d'entropie interne au système en plus de l'entropie d'échange, de sorte qu'en vertu de (2.7), cette entropie interne est non-négative :

$$\dot{S}_{\text{int}} \geq 0 \quad ; \quad \dot{S} = \dot{S}_{\text{ext}} + \dot{S}_{\text{int}} \qquad (2.8)$$

ce qui revient à écrire (2.7) comme :

$$\dot{S} \geq \frac{\dot{Q}}{T_0} \qquad (2.9)$$

Pour un milieu continu, le second principe de la thermodynamique s'exprime comme *l'inégalité fondamentale de l'entropie* :

$$\frac{d}{dt} \int_{\Omega} \rho s \, d\Omega \geq \int_{\Omega} \frac{r}{T} \, d\Omega - \int_{\partial\Omega} \frac{\mathbf{q} \cdot \mathbf{n}}{T} \, dS \qquad (2.10)$$

où s est l'entropie spécifique (par unité de masse). Le second principe fournit une borne supérieure du taux de quantité de chaleur que peut recevoir le volume Ω à une température T. L'inégalité (2.10) est à rapprocher de (2.9), considérant la température T uniforme et égale à T₀.

Le second principe s'écrit aussi sous forme locale, en faisant apparaître le taux de production irréversible d'entropie volumique $\dot{\Delta}$ :

$$\dot{\Delta} = \rho \dot{s} + \text{div} \left( \frac{\mathbf{q}}{T} \right) - \frac{r}{T} \geq 0 \qquad (2.11)$$

cette dernière devant donc être non-négative. L'inégalité (2.11) est la traduction locale de (2.8).

<footnote n="3">³ Il s'agit par exemple de la transformation d'une énergie mécanique en chaleur.</footnote>
## Inégalité de Clausius-Duhem

En combinant (2.6) et (2.11), il vient une nouvelle forme du second principe, connue sous le nom d'inégalité de Clausius-Duhem :

$$ \mathcal{D} = T \dot{\Delta} = \sigma : \mathbf{D} + \rho(T\dot{s} - \dot{e}) - \frac{\mathbf{q}}{T} \cdot \nabla T \ge 0 \quad (2.12) $$

où $\mathcal{D}$ est la dissipation volumique totale. Cette inégalité peut aussi s'écrire en fonction de l'énergie libre spécifique $\psi$, définie telle que :

$$ \psi = e - Ts \quad (2.13) $$

de sorte que l'inégalité de Clausius-Duhem se réécrive :

$$ \mathcal{D} = T \dot{\Delta} = \sigma : \mathbf{D} - \rho(\dot{\psi} + s\dot{T}) - \frac{\mathbf{q}}{T} \cdot \nabla T \ge 0 \quad (2.14) $$

La dissipation totale $\mathcal{D}$ se compose de deux termes correspondants à des phénomènes dissipatifs de nature physique différente, la première est mécanique (ou intrinsèque) $\mathcal{D}_{\text{int}}$, la seconde est thermique $\mathcal{D}_{\text{th}}$ :

$$ \mathcal{D}_{\text{int}} = \sigma : \mathbf{D} - \rho(\dot{\psi} + s\dot{T}) ; \quad \mathcal{D}_{\text{th}} = -\frac{\mathbf{q}}{T} \cdot \nabla T \quad (2.15) $$

On effectue généralement l'hypothèse de *découplage* de ces deux types de dissipation, de sorte qu'en assurant chacune d'être non-négative, leur somme soit aussi effectivement non-négative.

### 2.1.3 Écriture des lois de comportement

Le décompte lois de bilan/inconnues fait apparaître un surplus de onze inconnues, aussi il est nécessaire d'apporter onze équations additionnelles que sont les lois de comportement. Ces lois sont des modèles, plus des principes ou lois de bilan, elles traduisent la nature précise du milieu considéré. Les lois de bilan ne peuvent rien nous apprendre sur le comportement interne du milieu : c'est l'objet des lois de comportement que de rendre compte de la spécificité du milieu considéré. L'objet du cadre thermodynamique est d'offrir un cadre conceptuel pour l'écriture de lois de comportement qui, jointes aux équations de bilan, conduiront à un système fermé de lois d'évolution.

Bien qu'il y ait une grande latitudine pour construire ces lois, elles doivent cependant observer quelques restrictions. Les critères les plus communs sont la causalité (déterminisme), le principe d'objectivité (indépendance à l'observateur) et bien sûr elles doivent conduire à des processus thermodynamiquement admissibles. En effet, le second principe de la thermodynamique exprimant un principe d'évolution impose des restrictions à la construction des lois de comportement.

À ces restrictions, on ajoute pour une grande partie des modèles (mais pas tous !) le principe de l'état local et celui de la simplicité matérielle (théorie du premier gradient). La théorie de la plasticité qu'on développe ici entre dans ce cadre.

## 2.2 Méthode de l'état local

La méthode de l'état local postule que l'état thermomécanique d'un milieu matériel en un point et à un instant donnés est complètement défini par la connaissance des valeurs, à cet instant, d'un certain nombre de variables ne dépendant que du point considéré (c'est le principe de l'action locale). Les dérivées temporelles n'intervenant pas pour définir l'état, cette hypothèse
<figure><img src="image_35.png" /></figure>

implique que toute évolution puisse être considérée comme une succession d'états d'équilibre. Les phénomènes très rapides pour lesquels les échelles de temps seraient de l'ordre des temps de relaxation de retour à l'équilibre thermostatique sont donc exclus du champ d'application de cette théorie.

C'est par le choix de la nature et du nombre des variables d'état que l'on décrit plus ou moins finement les phénomènes physiques. Un processus sera qualifié de thermodynamiquement admissible si à chaque instant de l'évolution, le second principe est satisfait. Ces variables d'état sont généralement de deux sortes :

- les variables *observables* : la température T et la déformation totale $\boldsymbol{\epsilon}$, et intervenant dans les lois de bilan;

- les variables *internes*, associées aux phénomènes dissipatifs. L'état actuel dépend de l'histoire passée, représentée par les valeurs à chaque instant de ces variables internes. Dans le cas de la plasticité, on introduit le tenseur des déformations plastiques $\boldsymbol{\epsilon}^p$, mais aussi d'autres variables qui seront associées aux différents types d'écrouissage pris en compte. Ces variables seront désignées par $\mathbf{V}_1, \mathbf{V}_2, \dots, \mathbf{V}_k, \mathbf{V}_k$ pouvant représenter une variable scalaire ou tensorielle.

### 2.2.1 Potentiel thermodynamique et lois d'état

Une manière de décrire complètement l'état thermodynamique du milieu est de se donner un *potentiel thermodynamique*, à partir duquel se déduisent les *lois d'état*. L'énergie interne $e$ est un potentiel thermodynamique; par la transformation de Legendre, on peut en construire d'autres, en particulier l'énergie libre d'Helmholtz $\psi$ qu'on utilise couramment en mécanique des solides. L'énergie libre est une fonction des variables d'état (observables et internes); elle doit être concave par rapport à la température, et convexe par rapport aux autres variables d'état pour vérifier a priori les conditions de stabilité thermodynamiques.

$$ \psi = \psi(\boldsymbol{\epsilon}, T; \boldsymbol{\epsilon}^p, \mathbf{V}_k) \qquad (2.16) $$

où $\mathbf{V}_k$ représente la $k^e$ variable interne additionnelle considérée dans le modèle, $k = 1, \dots, n$. Les déformations n'intervenant en élastoplasticité que sous la forme de leur partition additive, (2.16) peut aussi s'écrire :

$$ \psi = \psi(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^p, T; \mathbf{V}_k) = \psi(\boldsymbol{\epsilon}^e, T; \mathbf{V}_k) \qquad (2.17) $$

soit, en différenciant le potentiel thermodynamique :

$$ \dot{\psi} = \frac{\partial \psi}{\partial \boldsymbol{\epsilon}^e} : \dot{\boldsymbol{\epsilon}}^e + \frac{\partial \psi}{\partial T} \dot{T} + \frac{\partial \psi}{\partial \mathbf{V}_k} \cdot \dot{\mathbf{V}}_k \qquad (2.18) $$

En reportant (2.18) dans l'expression de la dissipation intrinsèque (2.15), il vient :

$$ \left(\boldsymbol{\sigma} - \rho \frac{\partial \psi}{\partial \boldsymbol{\epsilon}^e}\right) : \dot{\boldsymbol{\epsilon}}^e + \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^p - \rho \left(s + \frac{\partial \psi}{\partial T}\right) \dot{T} - \frac{\partial \psi}{\partial \mathbf{V}_k} \cdot \dot{\mathbf{V}}_k \ge 0 \qquad (2.19) $$

En imaginant des transformations particulières, d'abord élastique à température constante :

$$ \dot{\boldsymbol{\epsilon}}^p = \dot{\mathbf{V}}_k = \dot{T} = 0, \quad \forall \boldsymbol{\epsilon}^e \qquad (2.20) $$

puis thermique de dilation n'engendrant aucune plasticité :

$$ \dot{\boldsymbol{\epsilon}}^p = \dot{\mathbf{V}}_k = 0, \quad \forall \dot{T} \qquad (2.21) $$
les lois d'état émergent naturellement :

$$
\begin{align}
\boldsymbol{\sigma} &= \rho \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^e} = \rho \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}} = -\rho \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^P} \tag{2.22} \\
s &= -\frac{\partial \psi}{\partial T} \nonumber
\end{align}
$$

Ces expressions définissent les lois de la thermoélasticité. Par analogie avec les relations précédentes, on définit les forces thermodynamiques $\mathbf{A}_k$ associées aux variables internes $\mathbf{V}_k$ par :

$$
\mathbf{A}_k = \rho \frac{\partial \psi}{\partial \mathbf{V}_k} \tag{2.23}
$$

Ce sont ces forces thermodynamiques qui gouvernent l'évolution des phénomènes irréversibles.
Tenant compte des lois d'état écrites, la dissipation intrinsèque se réécrit alors :

$$
D_{\text{int}} = \sigma : \dot{\boldsymbol{\epsilon}}^p - \mathbf{A}_k \cdot \dot{\mathbf{V}}_k \ge 0 \quad (2.24)
$$

L'inégalité (2.24) peut aussi se réécrire de façon à faire apparaître l'énergie mécanique fournie
au matériau lors du processus et l'énergie récupérable :

$$
D_{\text{int}} = \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} - (\boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^e + \mathbf{A}_k \cdot \dot{\mathbf{V}}_k) \ge 0 \quad (2.25)
$$

où $\boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}$ est la puissance mécanique fournie et $\boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^e + \mathbf{A}_k \cdot \dot{\mathbf{V}}_k$, la puissance récupérable. En résumé, les lois d'état sont construites en se donnant un potentiel thermodynamique. En général, un potentiel thermodynamique n'a pas de signification physique directe, on ne le mesure pas. Ce qu'on peut atteindre par l'expérience, ce sont les lois d'état. Les forces thermodynamiques forment un ensemble de variables normales en dualité avec les variables d'état, le tableau 2.1 liste ces variables introduites.

<table><caption>**Tableau 2.1** – Variables et forces thermodynamiques associées</caption>
  <thead>
    <tr>
      <th colspan="2">Variables d'état</th>
      <th rowspan="2">Forces thermodynamiques associées</th>
    </tr>
    <tr>
      <th>observables</th>
      <th>internes</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td><i>T</i></td>
      <td></td>
      <td><i>s</i></td>
    </tr>
    <tr>
      <td><i>&epsilon;</i></td>
      <td></td>
      <td><b>&sigma;</b></td>
    </tr>
    <tr>
      <td></td>
      <td><i>&epsilon;</i><sup>P</sup></td>
      <td>-<i>&sigma;</i></td>
    </tr>
    <tr>
      <td></td>
      <td><b>V</b><sub>k</sub></td>
      <td><b>A</b><sub>k</sub></td>
    </tr>
  </tbody>
</table>

2.2.2 Potentiel de dissipation et lois complémentaires

Les lois d'état sont insuffisantes pour décrire complètement le milieu. Il nous faut des lois
qui permettent de décrire l'avancement des phénomènes irréversibles en fonction des forces
thermodynamiques correspondantes : ce sont les *lois complémentaires* (ou d'évolution). La
seule restriction qu'elles peuvent prendre est d'assurer la positivité des dissipations intrinsèque
et thermique.

Pour définir ces lois relatives aux processus dissipatifs, on postule l'existence d'un potentiel,
ou pseudo-potentiel, de dissipation $\varphi$ s'exprimant comme une fonction à valeur scalaire continue
et convexe par rapport aux variables flux (et éventuellement des variables d'état), ou plutôt —
par commodité — de son dual $\varphi^*$, obtenu par la transformée de Legendre-Fenchel, de sorte
<figure><img src="image_36.png" /></figure>

que les lois complémentaires permettent d'obtenir les variables de type flux ($\dot{\mathbf{V}}_k$) en fonction des forces thermodynamiques (et aussi éventuellement des variables d'état) :

$$ \varphi(\dot{\mathbf{\epsilon}}^P, \dot{\mathbf{V}}_k, \mathbf{q}/T) \rightarrow \varphi^*(\sigma, \mathbf{A}_k, \nabla T) = \sup_{(\dot{\mathbf{\epsilon}}^P, \dot{\mathbf{V}}_k, \mathbf{q}/T)} (\mathcal{D} - \varphi(\dot{\mathbf{\epsilon}}^P, \dot{\mathbf{V}}_k, \mathbf{q}/T)) \quad (2.26) $$

où $\varphi^*$ doit être une fonction convexe, non-négative et nulle à l'origine :

$$ \varphi^* \ge 0 \quad ; \quad \varphi^*(0) = 0 \qquad (2.27) $$

Les lois complémentaires se déduisent alors de ce potentiel de dissipation en appliquant la propriété de normalité :

$$ \begin{aligned} \dot{\mathbf{\epsilon}}^P &= \frac{\partial \varphi^*}{\partial \sigma} \\ -\dot{\mathbf{V}}_k &= \frac{\partial \varphi^*}{\partial \mathbf{A}_k}, \quad k = 1, \dots, m \\ -\frac{\mathbf{q}}{T} &= \frac{\partial \varphi^*}{\partial (\nabla T)} \end{aligned} \qquad (2.28) $$

Nous dirons qu'un milieu est un *milieu standard généralisé* si il existe un (pseudo-)potentiel de dissipation $\varphi$ *convexe de ses arguments*, et minimum en zéro. La convexité de $\varphi$ assure le respect du second principe de la thermodynamique.

Pour un milieu standard généralisé, tout le problème de la modélisation des phénomènes réside dans la détermination de l'expression analytique du potentiel thermodynamique $\psi$ et du potentiel de dissipation $\varphi$ ou de son dual $\varphi^*$. La détermination de ces deux fonctions scalaires est nettement plus simple que l'identification de $m + 6 + 3$ lois complémentaires liant forces et vitesses généralisées. On parle aussi de thermodynamique des milieux à dissipation normale, au sens où on applique le principe de normalité au (pseudo-)potentiel de dissipation choisi pour déduire les lois complémentaires.

La construction d'une théorie de plasticité nécessite généralement la définition d'une fonction de charge (le critère), d'une fonction d'écoulement (qui donnera la direction de l'écoulement plastique en appliquant le principe de normalité), et d'une fonction d'écrouissage (même raisonnement pour l'écrouissage). Dans le cas d'un *matériau standard généralisé*, ces trois entités sont données par un seul potentiel : le (pseudo-)potentiel de dissipation $\varphi^*$.

### 2.2.3 Thermodynamique des processus irréversibles

Une première simplification consiste à supposer que le pseudo-potentiel de dissipation $\varphi^*$ est une forme quadratique définie positive des forces thermodynamiques, donnée par :

$$ \varphi^* = \frac{1}{2} \mathbf{C}_{\alpha\beta} \mathbf{A}_{\alpha} \mathbf{A}_{\beta} \qquad (2.29) $$

de sorte que les variables flux soient des combinaisons linéaires des forces thermodynamiques :

$$ \dot{\mathbf{V}}_{\alpha} = \mathbf{C}_{\alpha\beta} \mathbf{A}_{\beta} \qquad (2.30) $$

Les coefficients $\mathbf{C}_{\alpha\beta}$ sont les coefficients d'Onsager. L'hypothèse d'Onsager consiste à considérer la matrice $\mathbf{C}_{\alpha\beta}$ symétrique.

La thermodynamique des processus irréversibles (TPI) apparaît comme un cas particulier de la thermodynamique des processus à dissipation normale, la première étant une version linéaire de la seconde.
## 2.3 Plasticité associée

### 2.3.1 Principe de la maximisation de la dissipation intrinsèque

En plasticité, le comportement est indépendant de la vitesse de chargement. Aussi, il faut assurer qu'au cours d'une charge plastique, l'état de contrainte actuel reste sur la surface de charge $f = 0$, de sorte que l'état $f > 0$ soit impossible. Ceci impose quelques restrictions sur le potentiel $\phi^*$ ; en effet, ce potentiel définit d'une façon générale des surfaces équipotentielles dans l'espace des variables flux, associées à différentes vitesses de chargement. Ce cas général englobe donc la possible prise en compte d'effets visqueux dans le cadre des matériaux standards généralisés, notamment dans le cas de la viscoplasticité.

De ce point de vue, la plasticité apparaît comme un cas limite de la viscoplasticité, dans lequel ces surfaces équipotentielles vont venir se superposer à la surface d'écoulement, et donc assurer l'indépendance à la vitesse de chargement. Ce cas limite est obtenu en introduisant le principe de la maximisation de la dissipation intrinsèque sous contrainte que $f$ reste négatif ou nul :

$$ \mathcal{D}_{\text{int}}(\boldsymbol{\sigma}, \mathbf{A}_k) = \max_{(\boldsymbol{\sigma}^*, \mathbf{A}_k^*)} (\boldsymbol{\sigma}^* : \dot{\boldsymbol{\epsilon}}^P - \mathbf{A}_k^* \cdot \dot{\mathbf{V}}_k), \quad \forall (\boldsymbol{\sigma}^*, \mathbf{A}_k^*) \in (\mathbb{P} \times \mathcal{A}^{\text{ad}}) \quad (2.31) $$

La maximisation de la dissipation intrinsèque peut se transformer en un problème de minimisation de l'opposé de cette dissipation, et en introduisant un multiplicateur de Lagrange $\dot{\lambda}$ pour traiter la contrainte de type égalité, de sorte que :

$$ -\mathcal{D}_{\text{int}}(\boldsymbol{\sigma}, \mathbf{A}_k) = \min_{(\boldsymbol{\sigma}^*, \mathbf{A}_k^*)} (-\boldsymbol{\sigma}^* : \dot{\boldsymbol{\epsilon}}^P + \mathbf{A}_k^* \cdot \dot{\mathbf{V}}_k + \dot{\lambda} f), \quad \forall (\boldsymbol{\sigma}^*, \mathbf{A}_k^*) \in (S^{\text{ad}} \times \mathcal{A}^{\text{ad}}) \quad (2.32) $$

Cette minimisation conduit aux lois complémentaires de la plasticité associée :

$$ \dot{\boldsymbol{\epsilon}}^P = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} \quad (2.33) $$

$$ -\dot{\mathbf{V}}_k = \dot{\lambda} \frac{\partial f}{\partial \mathbf{A}_k} \quad (2.34) $$

Le principe de la maximisation de la dissipation intrinsèque a pour conséquences :

- l'introduction du multiplicateur plastique (qui est en fait un multiplicateur de Lagrange). Celui-ci étant indéterminé dans la loi de comportement, on introduit alors la condition de cohérence $\dot{f} = 0$, qui permettra sa détermination;

- la direction d'écoulement est normale à la surface de charge $f$ ; le domaine défini par $f$ est convexe. La définition de la direction d'écoulement à partir de la surface de charge conduit à la plasticité associée, au sens associée au critère de plasticité;

- l'existence des conditions de charge-décharge dites conditions complémentaires de Kühn-Tucker :

$$ \dot{\lambda} \ge 0 \ ; \ f \le 0 \ ; \ \dot{\lambda} f = 0 \quad (2.35) $$

### 2.3.2 Particularisation au critère de von Mises

Dans la suite, on utilisera le critère de von Mises<fnref n="4" /> dont la définition repose sur le second invariant du déviateur des contraintes. Pour un matériau élastique parfaitement plastique, ce critère s'écrit :

$$ f = \sigma_{\text{eq}} - \sigma_y \le 0 \quad \text{avec} \quad \sigma_{\text{eq}} = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} \quad \text{et} \quad \mathbf{s} = \boldsymbol{\sigma} - \frac{1}{3}(\text{tr} \, \boldsymbol{\sigma})\mathbf{1} \quad (2.36) $$

<footnote n="4">⁴ Aussi appelée la *J₂-flow theory* en anglais.</footnote>
<figure><img src="image_37.png" /></figure>

Comme évoqué au cours du premier chapitre, la régularité de la surface de charge définie par ce critère permet d'éviter des points singuliers pour lesquels la normale est définie dans un cône convexe, et donc des complications supplémentaires inutiles. L'application de la loi d'écoulement (2.33) à ce critère conduit à :

$$ \dot{\boldsymbol{\epsilon}}^P = \frac{3}{2} \dot{\lambda} \frac{\boldsymbol{s}}{\sigma_{eq}} \qquad (2.37) $$

On observe ainsi que pour un matériau régi par le critère de von Mises, la vitesse de déformation plastique est colinéaire au déviateur des contraintes. On vérifie aussi l'incompressibilité des déformations plastiques, le déviateur des contraintes étant à trace nulle par définition.

On notera d'autre part que la contrainte équivalente est définie avec un facteur 3/2 tandis que la vitesse de déformation plastique cumulée est définie avec un facteur 2/3. Ces définitions permettent de garder la cohérence avec la contrainte et la déformation dans le cas uniaxial, et de conserver l'équivalence énergétique entre l'écriture multiaxiale de la puissance plastique et l'écriture « équivalente » scalaire :

$$ \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^P = \sigma_{eq} \dot{\boldsymbol{p}} \qquad (2.38) $$

### 2.3.3 Matériaux élastique et parfaitement plastique

Le cas du matériau élastique parfaitement plastique est caractérisé par l'absence d'écrouissage. Ce matériau est décrit par le critère et le potentiel thermodynamique suivants :

$$ f = \sigma_{eq} - \sigma_y \le 0 \qquad (2.39) $$

$$ \rho\psi = \frac{1}{2} \mathbf{K} : (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^P) : (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^P) \qquad (2.40) $$

où $\mathbf{K}$ est le tenseur d'élasticité. On en déduit la loi de comportement qui associe la loi d'élasticité et la loi d'écoulement plastique :

$$ \boldsymbol{\sigma} = \mathbf{K} : (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^P) \\ \dot{\boldsymbol{\epsilon}}^P = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} \quad ; \quad \dot{\lambda} \ge 0 \quad ; \quad f \le 0 \quad ; \quad \dot{\lambda} f = 0 \qquad (2.41) $$

soit en combinant les deux :

$$ \dot{\boldsymbol{\epsilon}} = \mathbf{K}^{-1} : \dot{\boldsymbol{\sigma}} + \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} \qquad (2.42) $$

Le multiplicateur plastique est indéterminé pour ce type de matériau en vitesse de contrainte imposée. En effet, le module plastique étant nul (plateau dans la courbe contrainte-déformation), il existe une infinité de positions équivalentes en déformation plastique pour un état de contrainte admissible donné. Aussi, comme évoqué au premier chapitre, il est déterminé par la résolution du problème aux limites.

Le fait d'imposer la vitesse de déformation totale permet de lever cette indétermination, le multiplicateur plastique peut alors être déterminé en combinant la loi d'élasticité écrite en vitesse et la condition de cohérence :

$$ \dot{\boldsymbol{\sigma}} = \mathbf{K} : (\dot{\boldsymbol{\epsilon}} - \dot{\boldsymbol{\epsilon}}^P) \quad \text{et} \quad \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} = 0 \qquad (2.43) $$

Combinant les équations (2.43) et la loi d'écoulement plastique, il vient l'expression du multiplicateur plastique pour un matériau élastique parfaitement plastique<fnref n="5" /> :

$$ \dot{\lambda} = \left\langle \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbf{K} : \dot{\boldsymbol{\epsilon}} \quad ; \quad \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbf{K} : \frac{\partial f}{\partial \boldsymbol{\sigma}} \right\rangle^{+} \qquad (2.44) $$

<footnote n="5">⁵ La notation $\langle \cdot \rangle^+$ indique la partie positive de la quantité $\cdot$.</footnote>
Dans le cas particulier de l'élasticité isotrope et du critère de von Mises, l'expression du multiplicateur se simplifie et donne :

$$ \dot{\lambda} = \left\langle \frac{2}{3} \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\varepsilon}} \right\rangle^{+} \qquad (2.45) $$

ce qui donne pour un chargement uniaxial ($\dot{\varepsilon} = \dot{\varepsilon}_{11}$) :

$$ \dot{\lambda} = \dot{\varepsilon} \operatorname{sgn}(\sigma) \qquad (2.46) $$

autrement dit :

$$ \dot{\varepsilon}^{\mathrm{P}} = \dot{\varepsilon} \qquad (2.47) $$

Il est intéressant d'observer qu'on retrouve par le cadre thermodynamique la même loi de comportement élastoplastique que celle déduite du principe de Hill au premier chapitre, dans le cas d'un matériau élastique parfaitement plastique. En effet, en l'absence d'écrouissage, la maximisation de la dissipation intrinsèque se réduit à la maximisation de la dissipation plastique, c'est-à-dire en fait, le principe de Hill.

### 2.3.4 Lois à écrouissage isotrope

L'introduction de l'écrouissage dans le modèle se traduit par l'introduction d'une variable interne. Comme évoqué au premier chapitre, la variable introduite dans le cas de l'écrouissage isotrope doit permettre de paramétrer l'expansion du domaine élastique (figure 1.13). On utilise classiquement la déformation plastique cumulée $p$ définie par (1.21) comme variable interne en plasticité associée. Cette variable, introduite au sein du critère, permet de décrire la dilatation du domaine élastique :

$$ f(\boldsymbol{\sigma}, R) = \sigma_{\mathrm{eq}} - R(p) - \sigma_y \le 0 \qquad (2.48) $$

où $R$ est la force thermodynamique associée à la déformation plastique cumulée, définie par la loi d'état :

$$ R = \rho \frac{\partial \psi}{\partial p} \quad \text{avec} \quad R(p = 0) = 0 \qquad (2.49) $$

On choisit généralement un potentiel thermodynamique de la forme :

$$ \rho\psi(\boldsymbol{\epsilon}, p) = \rho\psi^e(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^P) + \rho\psi^P(p) \qquad (2.50) $$

où $\psi^e$ est l'énergie élastique emmagasinée donnée par (2.40), $\psi^P$ est l'énergie réversible stockée par modification de l'état interne du matériau$^6$. On observe que l'énergie libre ne se réduit pas en général à l'énergie élastique emmagasinée. Notons aussi que la déformation $\boldsymbol{\epsilon}$ n'intervient que dans l'énergie élastique emmagasinée $\psi^e$.

Le critère permet d'établir la loi d'évolution de la déformation plastique cumulée :

$$ \dot{p} = -\dot{\lambda} \frac{\partial f}{\partial R} = \dot{\lambda} = \sqrt{\frac{2}{3} \dot{\varepsilon}^P : \dot{\varepsilon}^P} \qquad (2.51) $$

Quelle que soit la forme choisie pour $R(p)$, la condition de cohérence permet de trouver le multiplicateur plastique, et s'écrit :

$$ \dot{f} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + \frac{\partial f}{\partial p} \dot{p} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} - \frac{\partial R}{\partial p} \dot{p} = 0 \qquad (2.52) $$

<footnote>$^6$ Il s'agit, par exemple, d'une énergie élastique stockée en raison de contraintes résiduelles microscopiques.</footnote>
<figure><img src="image_38.png" /></figure>

Si on suppose satisfaite l'hypothèse d'écrouissage positif ($\partial R / \partial p > 0$), il est aisé de déterminer l'expression du multiplicateur plastique :

$$ \dot{\lambda} = \left\langle \frac{\frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}}}{\frac{\partial R}{\partial p}} \right\rangle^{+} \qquad (2.53) $$

La loi obtenue en utilisant le critère de von Mises et une règle d'écrouissage isotrope est connue comme la loi de Prandtl-Reuss. Ce type de loi fonctionne bien essentiellement pour les chargements monotones.

### 2.3.5 Lois à écrouissage cinématique linéaire

Cet écrouissage est connu sous le nom de la loi de Prager. On introduit la variable d'écrouissage $\mathbf{X}$ (tenseur déviateur du second ordre), qui décrit la translation du centre du domaine d'élasticité (figure 1.14), dans le critère :

$$ f(\boldsymbol{\sigma}, \mathbf{X}) = \sigma_{\text{eq}}(\boldsymbol{\sigma} - \mathbf{X}) - \sigma_y \le 0 \quad \text{avec} \quad \sigma_{\text{eq}}(\boldsymbol{\sigma} - \mathbf{X}) = \sqrt{\frac{3}{2}(\mathbf{s} - \mathbf{X}) : (\mathbf{s} - \mathbf{X})} \quad (2.54) $$

Avec le critère de von Mises, la normale à la surface de charge s'écrit de la façon suivante :

$$ \mathbf{n} = \frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{3}{2} \frac{\mathbf{s} - \mathbf{X}}{\sigma_{\text{eq}}(\boldsymbol{\sigma} - \mathbf{X})} \qquad (2.55) $$

La variable d'écrouissage (ou force thermodynamique) $\mathbf{X}$ est associée à la variable interne $\boldsymbol{\alpha}$, de sorte qu'en définissant l'énergie libre telle que :

$$ \rho \psi(\boldsymbol{\varepsilon}, \boldsymbol{\alpha}) = \rho \psi^e(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) + \rho \psi^\alpha(\boldsymbol{\alpha}) \qquad (2.56) $$

avec :

$$ \rho \psi^e = \frac{1}{2} \mathbf{K} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) \quad ; \quad \rho \psi^\alpha(\boldsymbol{\alpha}) = \frac{1}{2} H \boldsymbol{\alpha} : \boldsymbol{\alpha} \qquad (2.57) $$

on en déduise les lois d'état et les lois complémentaires suivantes :

$$ \begin{array}{l} \boldsymbol{\sigma} = \mathbf{K} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) \\ \mathbf{X} = \rho \frac{\partial \psi}{\partial \boldsymbol{\alpha}} = H \boldsymbol{\alpha} \\ \left. \begin{array}{l} \dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} \\ \dot{\boldsymbol{\alpha}} = -\dot{\lambda} \frac{\partial f}{\partial \mathbf{X}} = \dot{\boldsymbol{\varepsilon}}^p \end{array} \right\} \quad \dot{\lambda} \ge 0 \ ; \quad f \le 0 \ ; \quad \dot{\lambda} f = 0 \end{array} \qquad (2.58) $$

De la même manière, la condition de cohérence permet la détermination du multiplicateur plastique, et s'écrit :

$$ \begin{align} \dot{f} &= \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + \frac{\partial f}{\partial \mathbf{X}} : \dot{\mathbf{X}} = 0 \notag \\ &= \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} - H \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\varepsilon}}^p \tag{2.59} \\ &= \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} - H \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} : \frac{\partial f}{\partial \boldsymbol{\sigma}} \notag \end{align} $$
soit dans l'hypothèse d'un module d'écrouissage $H$ positif, il vient :

$$ \dot{\lambda} = \left\langle \frac{\frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}}}{H \frac{\partial f}{\partial \boldsymbol{\sigma}} : \frac{\partial f}{\partial \boldsymbol{\sigma}}} \right\rangle^{+} \qquad (2.60) $$

L'écrouissage cinématique est généralement d'une grande importance en chargements cycliques, il est cependant difficile d'introduire une relation non linéaire entre $X$ et $\alpha$ dans le cadre présenté, puisque celle-ci conduirait inévitablement à une non-linéarité biunivoque non conforme avec les observations expérimentales. La figure 2.1 illustre à quoi conduit une relation non linéaire biunivoque pour l'écrouissage cinématique. Nous verrons à la section 2.4 comment lever ce

<figure><img src="image_39.png" /><figcaption>Figure 2.1 – Écrouissage cinématique : relations non linéaire (gauche) et linéaire (droite) entre $\alpha$ et $X$</figcaption></figure>

verrou.

En résumé, l'élastoplasticité associée avec normalité généralisée correspondant au matériau standard généralisé se définit à partir d'un potentiel thermodynamique (ici l'énergie libre spécifique $\psi$) et du critère de plasticité. Les lois d'état et les lois complémentaires constituant la loi de comportement du solide élastoplastique en sont alors déduites.

$$ \begin{aligned} \boldsymbol{\varepsilon} &= \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^p \\ \psi &= \psi(\boldsymbol{\varepsilon}^e, T; \mathbf{V}_k) \\ \boldsymbol{\sigma} &= \rho \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}} ; \quad s = -\frac{\partial \psi}{\partial T} ; \quad \mathbf{A}_k = \rho \frac{\partial \psi}{\partial \mathbf{V}_k} \\ f(\boldsymbol{\sigma}, \mathbf{A}_k) &\le 0 \\ \dot{\boldsymbol{\varepsilon}}^p &= \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} ; \quad -\dot{\mathbf{V}}_k = \dot{\lambda} \frac{\partial f}{\partial \mathbf{A}_k} \\ \dot{\lambda} &\ge 0 ; \quad \dot{\lambda} f = 0 ; \quad \dot{\lambda} \dot{f} = 0 \end{aligned} \qquad (2.61) $$

## 2.4 Plasticité non associée

L'hypothèse de normalité associée au critère de plasticité a été introduite d'abord en se fondant sur une observation expérimentale, c'est-à-dire la normalité de l'écoulement plastique pour les métaux, puis en *postulant* l'extension de cette normalité aux variables internes caractérisant les phénomènes dissipatifs. Ce cadre particulièrement attractif, ne permet cependant pas de pouvoir rendre compte correctement du comportement de certains matériaux et/ou pour certains types de sollicitations. Un exemple de situation où ce cadre atteint ses limites a trait à la modélisation
<figure><img src="image_40.png" /></figure>

l'écrouissage cinématique non linéaire, particulièrement important en chargement cyclique. Si la plasticité associée s'avère fonctionner relativement bien pour les chargements monotones, elle ne permet cependant pas de rendre compte correctement de l'écrouissage cinématique non linéaire en chargement cyclique. Il faut alors passer à un cadre dit « non associé ».

L'idée du cadre « non associé » est d'introduire un pseudo-potentiel $\mathcal{F}$ définissant une surface d'écoulement différente de la surface de charge définie par le critère de plasticité. Le principe de normalité généralisée est alors appliqué à cette surface d'écoulement.

$$ \dot{\mathbf{V}}_k = -\dot{\lambda} \frac{\partial \mathcal{F}}{\partial \mathbf{A}_k} ; \quad \dot{\lambda} \ge 0 \qquad (2.62) $$

La définition du matériau est alors réalisée au travers de trois fonctions scalaires : le potentiel thermodynamique $\psi$, le critère de plasticité $f$ et un pseudo-potentiel $\mathcal{F}$ définissant une surface d'écoulement. Ce pseudo-potentiel doit être défini de manière à respecter le second principe de la thermodynamique, c'est-à-dire :

— $\mathcal{F}$ doit être convexe par rapport aux forces thermodynamiques $\mathbf{A}_k$ et doit contenir l'origine ($\mathcal{F}(0) = 0$) :

$$ \mathcal{D}_{\text{int}} = -\mathbf{A}_k \cdot \dot{\mathbf{V}}_k = \dot{\lambda} \mathbf{A}_k \cdot \frac{\partial \mathcal{F}}{\partial \mathbf{A}_k} \ge 0 \qquad (2.63) $$

— $\mathcal{F}$ doit être positif dès que $\dot{\mathbf{V}}_k \neq 0$.

## 2.4.1 Calcul du multiplicateur plastique

Le calcul du multiplicateur plastique s'effectue à l'aide de l'équation de cohérence, s'écrivant de la façon suivante :

$$ 0 = \dot{f} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + \frac{\partial f}{\partial \mathbf{A}_k} \cdot \dot{\mathbf{A}}_k \qquad (2.64) $$

en développant la dérivée du terme de droite et en introduisant les lois d'évolutions (2.62), il vient :

$$ 0 = \dot{f} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} - \dot{\lambda} \frac{\partial f}{\partial \mathbf{A}_l} \cdot \frac{\partial \mathbf{A}_l}{\partial \mathbf{V}_k} \cdot \frac{\partial \mathcal{F}}{\partial \mathbf{A}_k} \qquad (2.65) $$

en introduisant la définition des forces thermodynamique (2.23) à partir du même potentiel, il vient l'expression du multiplicateur plastique :

$$ \dot{\lambda} = \left\langle \frac{\frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}}}{\frac{\partial f}{\partial \mathbf{A}_l} \cdot \frac{\partial^2 \rho \psi}{\partial \mathbf{V}_k \partial \mathbf{V}_l} \cdot \frac{\partial \mathcal{F}}{\partial \mathbf{A}_k}} \right\rangle^+ \qquad (2.66) $$

## 2.4.2 Particularisation

La particularisation d'un modèle de plasticité non associée s'effectue en choisissant le potentiel thermodynamique, le critère de plasticité et le pseudo-potentiel, par exemple de la façon suivante :

$$ \rho \psi(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^\mathrm{p}, T, \mathbf{V}) = \frac{1}{2}(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^\mathrm{p} - \boldsymbol{\epsilon}^\mathrm{th}) : \mathbf{K}(T) : (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^\mathrm{p} - \boldsymbol{\epsilon}^\mathrm{th}) + \rho \psi_r(r) + \rho \psi_\alpha(\boldsymbol{\alpha}) \\ f(\mathbf{A}_k; T) = \sigma_{\mathrm{eq}}(\boldsymbol{\sigma} - \mathbf{X}) - (R + \sigma_y) \qquad (2.67) $$

$$ \mathcal{F}(\mathbf{A}_k; T) = f(\mathbf{A}_k; T) + f_r(-R; T) + f_\alpha(-\mathbf{X}; T) $$
Pour un matériau élastoplastique à écrouissages isotrope et cinématique non linéaires, les formes précédentes peuvent être définies telle que :

$$ \rho\psi_r = \frac{1}{2} \frac{R_\infty}{b} r^2 \qquad f_r = \frac{1}{2} \frac{R^2}{R_\infty} \qquad (2.68) $$

$$ \rho\psi_\alpha = \frac{1}{2} C \alpha : \alpha \qquad f_\alpha = \frac{1}{2} \frac{\gamma}{C} \mathbf{X} : \mathbf{X} \qquad (2.69) $$

où $\gamma$, $R_\infty$, $b$ et $C$ sont des constantes matériau. Aussi, ce modèle tient compte de trois variables internes : la déformation plastique $\boldsymbol{\varepsilon}^P$, une variable scalaire $r$ associée à l'écrouissage isotrope, et une variable tensorielle $\boldsymbol{\alpha}$ associée à l'écrouissage cinématique. Le tableau 2.2 récapitule la dualité des variables internes et forces thermodynamiques introduites dans le modèle. Les lois

<table><caption>**Tableau 2.2** – Variables et forces thermodynamiques associées dans le modèle</caption><thead><tr><th colspan="2">Variables d'état</th><th rowspan="2">Forces thermodynamiques associées</th></tr><tr><th>observables</th><th>internes</th></tr></thead><tbody><tr><td><i>T</i></td><td></td><td><i>s</i></td></tr><tr><td><i>&epsilon;</i></td><td></td><td><b>&sigma;</b></td></tr><tr><td></td><td><i>&epsilon;</i><sup>P</sup></td><td>-<i>&sigma;</i></td></tr><tr><td></td><td><i>r</i></td><td><i>R</i></td></tr><tr><td></td><td><b>&alpha;</b></td><td><b>X</b></td></tr></tbody></table>

d'état se déduisent alors du potentiel thermodynamique :

$$ \boldsymbol{\sigma} = \mathbf{K} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^P - \boldsymbol{\varepsilon}^{\text{th}}) \qquad (2.70) $$

$$ R = \frac{R_{\infty}}{b} r \qquad (2.71) $$

$$ \mathbf{X} = C\boldsymbol{\alpha} \qquad (2.72) $$

Les lois d'évolution se déduisent quant à elles du pseudo-potentiel de dissipation, et donnent :

$$ \dot{\boldsymbol{\varepsilon}}^P = \frac{3}{2} \dot{\boldsymbol{\lambda}} \frac{\mathbf{s} - \mathbf{X}}{\sigma_{\text{eq}}(\boldsymbol{\sigma} - \mathbf{X})} \qquad (2.73) $$

$$ \dot{r} = \dot{\lambda} \left( 1 - \frac{R}{R_{\infty}} \right) \qquad (2.74) $$

$$ \boldsymbol{\alpha} = \frac{3}{2} \dot{\boldsymbol{\lambda}} \frac{\mathbf{s} - \mathbf{X}}{\sigma_{\text{eq}}(\boldsymbol{\sigma} - \mathbf{X})} - \dot{\boldsymbol{\lambda}} \frac{\gamma}{C} \mathbf{X} \qquad (2.75) $$

Notons aussi que, compte tenu de la spécification au critère de von Mises présenté à la section 2.3.2, il y a équivalence des écritures multiaxiale et « équivalente » de la puissance plastique, d'où l'égalité du multiplicateur plastique et de la vitesse de la déformation plastique cumulée :

$$ \text{Norme de Mises} \Rightarrow (\boldsymbol{\sigma} - \mathbf{X}) : \dot{\boldsymbol{\varepsilon}}^P = \sigma_{\text{eq}} \dot{p} \rightarrow \dot{p} = \dot{\lambda} \qquad (2.76) $$

En combinant d'une part (2.71) et (2.74), et d'autre part (2.72), (2.76) et (2.75), il vient les lois d'écrouissage isotrope et cinématique non linéaire :

$$ \dot{R} = \frac{1}{b}(R_{\infty} - R)\dot{p} \qquad (2.77) $$

$$ \dot{\mathbf{X}} = C\dot{\boldsymbol{\varepsilon}}^P - \gamma\mathbf{X}\dot{p} \qquad (2.78) $$
<figure><img src="image_41.png" /></figure>

La relation (2.78) est la loi d'Armstrong-Frédéricks, introduite en 1966, bien avant d'ailleurs qu'elle ne soit intégrée au cadre thermodynamique. On observe que l'inconvénient de la loi de Prager évoqué à la section 2.3.5 est levé en introduisant un terme de rappel traduisant un effet de mémoire évanescente du trajet de déformation (dépendance à $p$). De cette manière, on élimine la linéarité et la relation biunivoque entre l'écrouissage et la déformation plastique. La figure 2.2 montre l'évolution de ces écrouissages dans un cas unidimensionnel en fonction de la déformation plastique. Cette loi permet notamment de décrire l'effet « rochet ».

<figure><img src="image_42.png" /><figcaption>**Figure 2.2** – Évolution des écrouissages isotrope (à gauche) et cinématique (à droite) en fonction de la déformation plastique dans un cas unidimensionnel</figcaption></figure>
# Calcul des structures élastoplastiques

## 3.1 Un problème d'évolution

La détermination de l'état d'équilibre d'une structure sous un chargement donné s'effectue par la résolution du problème aux limites, c'est-à-dire par la résolution d'un système d'équations aux dérivées partielles (les équations d'équilibre) auquelles sont ajoutées les équations de compatibilité géométrique du champ de déformation, la loi de comportement et les conditions aux limites.

Dans le cas d'une structure composée d'un matériau élastoplastique, la réponse du système dépend du trajet de chargement du fait de l'irréversibilité du comportement, et pas seulement du chargement instantané comme c'est le cas en élasticité. La résolution du problème aux limites pour une structure élastoplastique se formule alors comme un problème d'évolution non linéaire. Cette évolution implique l'existence d'un état initial qu'il faut aussi préciser dans l'écriture du problème. La non-linéarité introduite dans le problème provient ici du comportement du matériau. D'autres types de non-linéarités peuvent intervenir dans un problème aux limites, notamment d'origine géométrique (grandes transformations) ou au travers des conditions aux limites (contact unilatéral, rayonnement, etc.).

La résolution de ce type de problème étant généralement complexe, on limite le cadre d'étude en faisant les hypothèses de travail suivantes :

- on se place dans le cadre des petites transformations.

- les matériaux étudiés ont un comportement élastoplastique parfaitement plastique ou à écrouissage positif.

- le chargement est imposé de façon suffisamment lente pour se placer dans l'approximation des transformations quasi-statiques et isothermes.

Ce chapitre est organisé en quatre sections. Dans la section 3.2, on écrit la formulation forte du problème élastoplastique. La section 3.3 résume les principaux résultats d'existence et d'unicité des problèmes d'évolution élastoplastiques. On présente ensuite quelques exemples de résolution de problèmes d'élastoplasticité en petites transformations dans la section 3.4, au travers de l'étude de la plasticité dans les barres et dans les poutres. Enfin, dans la section 3.5, on montre que l'irréversibilité du comportement élastoplastique génère des contraintes résiduelles, lesquelles dérivent de l'incompatibilité des déformations plastiques. Il existe alors aussi des déformations résiduelles compatibles, qui sont obtenues en sommant les déformations plastiques aux déformations élastiques générées par les contraintes résiduelles. On parle enfin de l'origine de l'écrouissage de structure, généralement observé sur la courbe de charge d'une structure.
## 3.2 Formulation forte du problème

On considère une structure occupant un domaine $\Omega$, de frontière régulière $\partial\Omega$ se décomposant en deux sous parties non chevauchantes (indépendantes du temps pour simplifier) :

$$ 
\begin{align} 
\partial\Omega &= \partial\Omega_u \cup \partial\Omega_F \tag{3.1} \\ 
\emptyset &= \partial\Omega_u \cap \partial\Omega_F \nonumber
\end{align} 
$$

où $\partial\Omega_u$ et $\partial\Omega_F$ réfèrent respectivement aux parties du bord sur lesquelles le déplacement et le vecteur contrainte sont imposés. Ce corps est composé d'un matériau élastoplastique standard (plasticité associée) à écrouissage nul ou positif.

**État initial** On se donne l'état initial du système $(\mathbf{u}_0, \boldsymbol{\sigma}_0, \boldsymbol{\varepsilon}_0^p : (\mathbf{A}_k^0, \mathbf{V}_k^0))$. Ces champs ne sont pas indépendants ; l'état de contrainte initial $\boldsymbol{\sigma}_0$ et le champ de déformations plastiques initiales doivent être compatibles entre eux, de sorte qu'il existe un champ de déplacement $\mathbf{u}_0$ cinématiquement admissible avec les données initiales en déplacement tel que :

$$ \boldsymbol{\varepsilon}(\mathbf{u}_0) = \mathbf{K}^{-1} : \boldsymbol{\sigma}_0 + \boldsymbol{\varepsilon}_0^p \quad (3.2) $$

D'autre part, l'état interne du matériau caractérisé par la suite variables internes $(\mathbf{V}_k)$ doit aussi être compatible avec les forces thermodynamiques associées au travers des lois d'état :

$$ \mathbf{A}_k^0 = \rho \frac{\partial \psi}{\partial \mathbf{V}_k^0} \quad (3.3) $$

L'état initial non nul pour les différents champs mécaniques dépend de l'histoire antérieure au cours de laquelle le corps a pu subir des sollicitations :

- cet état peut provenir d'un précédent trajet de chargement subi par la structure. Ceci est notamment fréquent en simulation numérique des procédés de fabrication et/ou de mise en forme comme le formage incrémental par exemple. La structure va passer par plusieurs phases de mise en forme, chacune déformant et modifiant l'état résiduel de la structure. Il règne alors dans celle-ci un champ d'autocontrainte, qui par définition, est compatible avec un chargement nul. Soit en appliquant le principe des puissance virtuelles, ce champ vérifie :

$$ \int_{\Omega} \boldsymbol{\sigma}_0 : \boldsymbol{\varepsilon} \, d\Omega = 0, \quad \forall \boldsymbol{\varepsilon} \text{ cinématiquement admissible} \quad (3.4) $$

- cet état peut aussi venir d'une précharge du système, auquel cas l'état initial $(\boldsymbol{\sigma}_0, \mathbf{u}_0)$ est obtenu en résolvant un problème d'élasticité linéaire avec précontraintes sous le chargement initial $(\mathbf{u}_0^d, \mathbf{F}_0^d, \mathbf{f}_0^d)$, de sorte que le problème s'écrive :

$$ 
\begin{align} 
\mathrm{div} \, \boldsymbol{\sigma}_0 + \mathbf{f}_0^d &= \mathbf{0} && \forall \mathbf{x} \in \Omega \nonumber \\ 
\boldsymbol{\sigma}_0 &= \mathbf{K} : (\boldsymbol{\varepsilon}_0 - \boldsymbol{\varepsilon}_0^p) && \forall \mathbf{x} \in \Omega \nonumber \\ 
\boldsymbol{\varepsilon}_0 &= (\nabla \mathbf{u}_0 + \nabla^\top \mathbf{u}_0)/2 && \forall \mathbf{x} \in \Omega \tag{3.5} \\ 
\mathbf{u}_0 &= \mathbf{u}_0^d && \forall \mathbf{x} \in \partial\Omega_u \nonumber \\ 
\boldsymbol{\sigma}_0 \cdot \mathbf{n} &= \mathbf{F}_0^d && \forall \mathbf{x} \in \partial\Omega_F \nonumber 
\end{align} 
$$

Le terme $-\mathbf{K} : \boldsymbol{\varepsilon}_0^p$ dans la loi d'élasticité joue le rôle d'une précontrainte. Il reste ensuite à vérifier le critère de plasticité avec cet état initial. On observe qu'un champ dit d'autocontrainte apparaît alors comme un cas particulier de celui-ci en fixant le chargement initial à zéro.
**Chargement** On se donne aussi l'histoire du chargement $\mathbf{f}^d(t)$, $\mathbf{F}^d(t)$, $\mathbf{u}^d(t)$ auquel est soumis le corps dans l'intervalle de temps $[0, T]$, dont on suppose l'évolution suffisamment régulière au regard des hypothèses précédemment énoncées. Du fait de l'indépendance du comportement vis-à-vis des vitesses, le temps est simplement ici un paramètre servant à décrire la séquence des états du chargement et de la structure.

**Problème aux limites** La résolution du problème aux limites consiste à trouver à tout instant $t \in [0, T]$ un champ de déplacement et un champ de contrainte (le couple $(\boldsymbol{\sigma}, \mathbf{u})$) respectivement cinématiquement et statiquement admissibles, et vérifiant la relation de comportement au travers des équations de compatibilité.

<figure><img src="image_43.png" /><figcaption>**Figure 3.1 – Problème aux limites**</figcaption></figure>

Le problème s'énonce alors : trouver $(\mathbf{u}, \boldsymbol{\sigma}, \boldsymbol{\varepsilon}^P; (\mathbf{A}_k, \mathbf{V}_k))$ sur $\Omega \times [0, T]$ tels que :

□ **u** cinématiquement admissible (C.A.)

* **u** « régulier » dans $\Omega$

* Conditions aux limites : $\mathbf{u} = \mathbf{u}^d, \forall \mathbf{x} \in \partial\Omega_u$

* Conditions initiales : $\mathbf{u}(\mathbf{x}, t=0) = \mathbf{u}_0(\mathbf{x})$

□ **σ** statiquement admissible (S.A.)

* Équilibre : $\text{div} \boldsymbol{\sigma} + \mathbf{f}^d = \mathbf{0}$ dans $\Omega$

* Conditions aux limites : $\mathbf{F} = \boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{F}^d, \forall \mathbf{x} \in \partial\Omega_F$

* Conditions initiales : $\boldsymbol{\sigma}(\mathbf{x}, 0) = \boldsymbol{\sigma}_0(\mathbf{x})$

□ **Compatibilité** (P)
$$ \boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla^\top \mathbf{u}) \quad \text{dans } \Omega $$

□ **Relation de comportement élastoplastique**

* Partition : $\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^P$

* Lois d'état : $\psi = \psi(\boldsymbol{\varepsilon}^e, T; \mathbf{V}_k)$, $\boldsymbol{\sigma} = \rho\partial_\boldsymbol{\varepsilon}\psi$, $s = -\partial_T\psi$, $\mathbf{A}_k = \rho\partial\mathbf{v}_k\psi$

* Lois d'évolution : $\dot{\boldsymbol{\varepsilon}}^P = \dot{\lambda}\partial_\boldsymbol{\sigma} f, -\dot{\mathbf{V}}_k = \dot{\lambda}\partial_{\mathbf{A}_k} f$

* Conditions Kühn-Tucker : $\dot{\lambda} \ge 0, f(\boldsymbol{\sigma}, \mathbf{A}_k) \le 0, \dot{\lambda} f = 0$

* Condition de cohérence : $\dot{\lambda} \dot{f} = 0$

* Conditions initiales : $\boldsymbol{\varepsilon}^P(\mathbf{x}, 0) = \boldsymbol{\varepsilon}_0^P(\mathbf{x}), \mathbf{A}_k(\mathbf{x}, 0) = \mathbf{A}_k^0(\mathbf{x}), \mathbf{V}_k(\mathbf{x}, 0) = \mathbf{V}_k^0(\mathbf{x})$

## 3.3 Résultats sur l'évolution des systèmes élastoplastiques

On peut se poser la question de l'existence et de l'unicité de la solution du problème d'évolution (P). Les résultats présentés sont valides dans le cadre de l'hypothèse des petites perturba-
tions. Il est montré[1] qu'en se donnant les hypothèses suivantes :

— le problème élastique linéaire posé avec les mêmes données que le problème élastoplastique admet une solution $(\mathbf{u}^{\text{el}}, \boldsymbol{\sigma}^{\text{el}})$ unique, à un éventuel mouvement de corps rigide près. Cela impose de la régularité sur le chargement imposé ainsi que la définie-positivité du tenseur d'élasticité $\mathbf{K}$ (existence d'une énergie interne $e$) et de son inverse le tenseur des souplesses $\mathbf{S}$.

— le matériau considéré est élastoplastique à écrouissage strictement positif.

Alors le problème d'évolution quasi-statique d'évolution d'une structure élastoplastique à écrouissage strictement positif admet une solution $(\mathbf{u}(\mathbf{x}, t), \boldsymbol{\sigma}(\mathbf{x}, t))$, $\forall t \in [0; T]$. Le champ de contrainte est unique, le champ de déplacement est donné à un mouvement de corps rigide près.

### Cas de la plasticité parfaite
La matériau élastique parfaitement plastique limite en chaque point du sytème la valeur du champ de contrainte par le critère de plasticité. Le champ de contrainte $\boldsymbol{\sigma}$ ne peut alors exister que si il est statiquement et plastiquement admissible. Sous cette condition, il est possible de montrer l'unicité du champ de contraintes pour ce problème[1,2,7]. Cependant, l'unicité du champ de déplacement n'est pas assurée à cause de la non-unicité du taux de déformation plastique $\dot{\boldsymbol{\sigma}}^P$. En combinant la partition du taux de déformation, la loi d'élasticité et la loi d'écoulement plastique, on peut écrire :

$$ \dot{\boldsymbol{\epsilon}} = \boldsymbol{S} : \dot{\boldsymbol{\sigma}} + \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} \qquad (3.6) $$

le multiplicateur plastique étant indéterminé (il est seulement astreint à être positif ou nul) dans le cas de la plasticité parfaite. La non-unicité du champ de déplacement traduit le phénomène d'écoulement plastique non contrôlable de la structure.

## 3.4 Problèmes d'élastoplasticité en petites transformations

On traite dans cette section deux problèmes d'élastoplasticité. On étudie dans un premier temps la plasticité dans les barres au travers d'une structure treillis composée de trois barres ; le comportement d'une barre étant unidimensionnel, ceci permet de simplifier la mise en œuvre. Dans un second temps, on étudie la plasticité dans les poutres avec le problème classique de la flexion trois points. Des problèmes tridimensionnels seront aussi abordés au cours des travaux dirigés.

### 3.4.1 Structure treillis à trois barres

On considère la structure treillis de la figure B.2 composée de trois barres, articulées entre elles au point d'application du chargement $F$, et chacune étant articulée par rapport au bâti. On

<figure><img src="image_44.png" /><figcaption>Figure 3.2 – Structure treillis de trois barres</figcaption></figure>
suppose que les barres sont à l'état initial sans contraintes, et que le chargement est appliqué suffisamment lentement pour rentrer dans le cadre établi précédemment. On suppose d'autre part que les trois barres ont un comportement élastique parfaitement plastique *identique*, c'est-à-dire un module d'Young E et une limite d'élasticité $\sigma_0$, et une section *identique* S. Pour simplifier l'étude, l'angle $\alpha$ sur la figure B.2 est égal à 45°.

**Phase élastique** Le problème est constitué de trois inconnues $N_1$, $N_2$ et $N_3$ représentant les efforts dans les barres, et de deux équations d'équilibre obtenues en isolant le nœud où la charge est appliquée (figure 3.3). Par conséquent, le système est hyperstatique de degré un. Les

<figure><img src="image_45.png" /><figcaption>**Figure 3.3 – Équilibre des forces**</figcaption></figure>

équations d'équilibre donnent :

$$N_1 = N_3 \quad ; \quad N_2 = F - \sqrt{2}N_1 \qquad (3.7)$$

L'inconnue hyperstatique peut être levée en considérant le comportement élastique des barres. Prenant $N_1$ comme inconnue hyperstatique, la résolution peut s'effectuer par exemple par une méthode énergétique, l'énergie de déformation élastique s'écrivant ici :

$$E_d = \frac{1}{2ES} (h\sqrt{2}N_1^2 + hN_2^2 + h\sqrt{2}N_3^2) \qquad (3.8)$$

soit en y introduisant les équations d'équilibre (3.7) :

$$E_d = \frac{h}{2ES} (2\sqrt{2}N_1^2 + (F - \sqrt{2}N_1)^2) \qquad (3.9)$$

Le théorème de Ménabréa implique alors :

$$\frac{\partial E_d}{\partial N_1} = 0 \Rightarrow N_1 = N_3 = \frac{F}{2 + \sqrt{2}} \quad \text{et} \quad N_2 = \frac{F\sqrt{2}}{1 + \sqrt{2}} \qquad (3.10)$$

**Domaine d'élasticité** Le domaine d'élasticité de la structure est engendré par celui des trois barres, de sorte que :

$$\{N_i \in \mathbb{R}^3, i = 1, 2, 3; f(N_i) = |N_i| - \sigma_0 S \le 0\} \qquad (3.11)$$

Ce domaine d'élasticité peut être représenté dans le plan $(N_1, N_2)$ compte tenu du fait que $N_1 = N_3$, par un carré de côté $2\sigma_0 S$ (figure 3.4). La solution élastique (3.10) peut être tracée dans le plan $(N_1, N_2)$ en éliminant F des équations (3.10), cette solution est alors représentée par la droite $N_2 = 2N_1$. La barre 2 sera donc la première à plastifier; la charge de première plastification de la structure, notée $F^{(1)}$, est alors donnée par :

$$F^{(1)} = \frac{1 + \sqrt{2}}{\sqrt{2}} \sigma_0 S \qquad (3.12)$$
<figure><img src="image_46.png" /><figcaption>Figure 3.4 – Domaine d’élasticité</figcaption></figure>

L'allongement de chacune des barres peut être calculé en fin de phase élastique avec la loi de comportement élastique d'une barre :

$$N_i = ES \frac{\delta_i}{L_i}, \quad i = 1, 2, 3 \qquad (3.13)$$

où $\delta_i$ dénote l'allongement de la $i$-e barre et $L_i$ sa longueur. L'allongement de la barre 2 est égal à celui du point d'application de la charge sur la structure, et est noté $\delta$ :

$$\delta = \delta_2 = \frac{h\sigma_0}{E} \qquad (3.14)$$

L'allongement des barres 1 et 3 est obtenu en utilisant la relation de compatibilité géométrique exprimant que les barres sont articulées au point d'application de la charge, soit $\varepsilon_2 = 2\varepsilon_1 = 2\varepsilon_3$ (on obtient la même chose par linéarité de (3.13)). Il vient l'allongement des barres 1 et 3 :

$$\delta_1 = \delta_3 = \frac{\delta}{\sqrt{2}} = \frac{h\sigma_0}{\sqrt{2}E} \qquad (3.15)$$

**Phase élastoplastique** On continue à présent de charger la structure au delà de $F^{(1)}$. Le matériau des barres étant élastique parfaitement plastique et la barre 2 étant plastifiée, on a toujours $N_2 = \sigma_0 S$. Les deux équations d'équilibre (3.7) permettent alors de déterminer les deux inconnues $N_1$ et $N_3$, soit :

$$N_3 = N_1 = \frac{F - \sigma_0 S}{\sqrt{2}} \qquad (3.16)$$

Il reste à vérifier si cette solution vérifie le critère de plasticité $f(N_i) = |N_i| - \sigma_0 S \le 0$ et la condition de charge $\Delta\varepsilon_2^p > 0$. Le critère appliqué aux barres 1 ou 3 nous donne la seconde charge de plastification $F^{(2)} = (1 + \sqrt{2})\sigma_0 S$, c'est-à-dire la charge pour laquelle les barres 1 et 3 plastifient. La condition de charge $\Delta\varepsilon_2^p > 0$ peut être vérifiée en utilisant la relation de compatibilité des déformations $\varepsilon_2 = 2\varepsilon_1$ et la loi de comportement des barres 1 et 3 qui restent élastiques durant cette phase :

$$\begin{gathered} \varepsilon_2 = 2\varepsilon_1 = \frac{2N_1}{ES} = \frac{\sqrt{2}}{ES}(F - \sigma_0 S) \rightarrow \delta_2 = \frac{\sqrt{2}h}{ES}(F - \sigma_0 S) \qquad (3.17) \\ \delta_1 = \delta_3 = \frac{h(F - \sigma_0 S)}{ES} \end{gathered}$$
En utilisant la partition de la déformation $\varepsilon_2 = \varepsilon_2^e + \varepsilon_2^p$, il vient :

$$ \varepsilon_2^p = \frac{1}{E} \left( \frac{F\sqrt{2}}{S} - (1 + \sqrt{2})\sigma_0 \right) > 0 \qquad (3.18) $$

On vérifie bien qu'à un incrément de charge $\Delta F > 0$ correspond un incrément positif de la déformation plastique de la barre 2, $\Delta\varepsilon_2^p > 0$. Cette solution est valable jusqu'à la charge de seconde plastification $F^{(2)}$. Au delà, toutes les barres sont plastiques, la charge ne peut plus augmenter (comportement élastique parfaitement plastique des barres), on dit qu'on a atteint la *charge limite* $F^{(L)} = F^{(2)}$. Une charge limite correspond à un état limite du système dans lequel il y a *écoulement plastique libre* (non limité) à charge constante. Il s'ensuit donc la *ruine de la structure*. Cette ruine est obtenue pour la flèche suivante de la structure :

$$ \delta = 2h\varepsilon_1 = \frac{2hN_1}{ES} = \frac{2\sigma_0h}{E} \qquad (3.19) $$

La figure 3.5 présente la courbe de charge de la structure. Il est intéressant d'observer qu'entre la charge de première plastification et la charge limite un « écrouissage de structure », bien que le matériau constitutif des barres soit élastique parfaitement plastique. Cet écrouissage est dû ici à la présence de zones encore non plastifiées, autrement dit c'est la présence de zones élastiques restantes qui génère cet effet de consolidation de la structure.

<figure><img src="image_47.png" /><figcaption>Figure 3.5 – Courbe de charge de la structure treillis</figcaption></figure>

**Décharge à partir de la charge limite** On note $\Delta F$ la variation de la charge $F$ à partir de la charge limite. Dans un premier temps, la décharge est supposée élastique. Aussi la variation des efforts dans les barres est donnée par la solution élastique (3.10), soit :

$$ \begin{aligned} \Delta N_1 &= \Delta N_3 = \frac{\Delta F}{2 + \sqrt{2}} \\ \Delta N_2 &= \frac{\Delta F \sqrt{2}}{1 + \sqrt{2}} \end{aligned} \qquad (3.20) $$

Or $\Delta F$ varie de $F^{(L)}$ à 0, donc $\Delta F = -F^{(L)} = -(1 + \sqrt{2})\sigma_0 S$, de sorte que les variations d'effort dans les barres se réécrivent :

$$ \begin{aligned} \Delta N_1 &= \Delta N_3 = -\frac{(1 + \sqrt{2})\sigma_0 S}{2 + \sqrt{2}} \\ \Delta N_2 &= -\sqrt{2}\sigma_0 S \end{aligned} \qquad (3.21) $$

La valeur des efforts dans les barres à la fin de la décharge est obtenue en sommant leurs valeurs à l'état limite ($N_1^{(L)} = N_2^{(L)} = N_3^{(L)} = \sigma_0 S$) et celles associées à la décharge élastique, il
vient :

$$N_1 = N_1^{(L)} + \Delta N_1 = N_3^{(L)} + \Delta N_3 = \frac{\sigma_0 S}{2 + \sqrt{2}} \quad (3.22)$$

$$N_2 = N_2^{(L)} + \Delta N_2 = (1 - \sqrt{2})\sigma_0 S$$

On observe qu'à la fin de la décharge il subsiste des efforts à l'intérieurs les barres ! Il existe donc des contraintes résiduelles après décharge. Ici, c'est l'hyperstaticité de la structure qui génère les contraintes résiduelles. Il est possible de vérifier que ces efforts résiduels satisfont les équations d'équilibre (3.7) puisque $N_2 = -\sqrt{2}N_1$. Aussi les contraintes résiduelles sont un champ autoéquilibré, c'est-à-dire statiquement admissible avec une charge nulle ($F = 0$). Les contraintes résiduelles proviennent du fait que les déformations plastiques ne vérifient pas les équations de compatibilité, en effet $\epsilon_2^P \neq 0$ et $\epsilon_1^P = \epsilon_3^P = 0$ entraînent que $\epsilon_2^P \neq 2\epsilon_1^P$.

À décharge, les barres 1 et 3 sont en traction, tandis que la barre 2 est passée en compression. Il reste à vérifier que la décharge est bien élastique, en vérifiant le critère pour chacune des barres :

$$\forall i, |N_i| - \sigma_0 S \le 0 \quad \text{avec } N_i = N_i|_{F=F(L)} + \Delta N_i \quad (3.23)$$

On peut calculer que c'est la barre 2 qui plastifiera la première en compression pour un accroissement de $\Delta N_2 = -2\sigma_0 S$ correspondant à une variation de charge de :

$$\Delta F = -(2 + \sqrt{2})\sigma_0 S \quad (3.24)$$

c'est-à-dire pour une variation de charge plus importante que celle opérée par la décharge ; la décharge est donc bien élastique. Enfin, la figure 3.6 présente l'évolution des efforts dans

<figure><img src="image_48.png" /><figcaption>Figure 3.6 – Efforts dans les barres</figcaption></figure>

chaque barre en fonction de la charge appliquée à la structure au cours du cycle charge-décharge. L'effort dans la barre 2 croît deux fois plus vite que dans les autres barres durant la phase élastique, en raison de la condition de compatibilité ($\epsilon_2 = 2\epsilon_1$) et de la linéarité de la loi élastique. Durant la phase élastoplastique, les efforts dans les barres 1 et 3 croissent plus vite que dans la phase élastique, l'état limite de la structure est atteint lorsqu'elles plastifient. Enfin, la condition de compatibilité et la décharge étant élastique impliquent que la barre 2 passe en compression, alors que les barres 1 et 3 restent en traction à décharge.

### 3.4.2 Exercices sur les barres en élastoplasticité

Résoudre le même problème :
1. Pour un chargement cyclique de la structure, de sorte que les état limites soient atteints en traction et en compression, toujours avec un matériau élastique parfaitement plastique.

2. Prendre en compte l'écrouissage dans les barres, d'abord cinématique puis isotrope, et effectuer le même cycle de chargement. Selon le type d'écrouissage pris en compte, le domaine d'élasticité de la structure variera comme présenté sur la figure 3.7.

<figure><img src="image_49.png" /><figcaption>Figure 3.7 – Modification du domaine d’élasticité de la structure treillis par écrouissage</figcaption></figure>

### 3.4.3 Plasticité dans les poutres en flexion simple

On considère le problème de flexion trois points présenté sur la figure 3.8. La poutre de section

<figure><img src="image_50.png" /><figcaption>Figure 3.8 – Flexion trois points</figcaption></figure>

symétrique est chargée en son milieu par une force supposée ponctuelle $F$ et repose sur deux appuis à chacune de ses extrémités. On suppose que la poutre obéit à l'hypothèse d'Euler-Bernoulli, qu'elle est à l'état initial sans contraintes et est constituée d'un matériau élastique parfaitement plastique, de module d'Young $E$ et de limite d'élasticité $\sigma_0$. On note $I$ le moment quadratique de la poutre par rapport à la direction $z$. Pour simplifier le problème, on considère que la poutre a une section rectangulaire de hauteur $h$ et de largeur $b$.

On note $\delta$ la flèche au milieu de la poutre. La solution élastique est connue et cette flèche est donnée par :

$$ \delta = \frac{FL^3}{48EI} \qquad (3.25) $$

Le torseur de cohésion (représentant les efforts intérieurs généralisés dans la poutre) est composé d'un effort tranchant dans la direction $y$ et d'un moment fléchissant autour de la direction $z$, dont le diagramme est présenté sur la figure 3.9. Compte tenu de la répartition de la contrainte normale dans la section pour une poutre soumise à de la flexion, les fibres les plus éloignées
<figure><img src="image_51.png" /><figcaption>**Figure 3.9** – Diagramme du moment fléchissant</figcaption></figure>

de la fibre neutre situées à une ordonnée $\pm h/2$ plastifieront en premier. En négligeant la contribution due au cisaillement, le critère de plasticité s'écrit :

$$f = |\sigma| - \sigma_0 \le 0; \quad \sigma = -\frac{M_{fz}}{I} y \qquad (3.26)$$

où y est l'ordonnée courante dans la section. En appliquant le critère, la charge de première plastification est donnée par :

$$F^{(1)} = \frac{2bh^2\sigma_0}{3L} \qquad (3.27)$$

soit une flèche de :

$$\delta^{(1)} = \frac{\sigma_0 L^2}{6Eh} \qquad (3.28)$$

**Phase élastoplastique** Les fibres extrémales de la poutre ayant plastifiées, la zone plastique va s'étendre vers la fibre moyenne de la poutre à mesure que la charge imposée va croître. La répartition de la contrainte normale dans la section n'est plus linéaire, mais atteint un palier à une ordonnée notée c sur la figure 3.10, dû au comportement parfaitement plastique de la poutre; en effet d'une façon générale, la plasticité limite les contraintes dans une structure. Avec cette répartition de contrainte dans la section plastifiée, le moment fléchissant se calcule

<figure><img src="image_52.png" /><figcaption>**Figure 3.10** – Contrainte normale dans la section en phase élastoplastique</figcaption></figure>

par définition :

$$ \begin{align} M_{fz} &= b \int_{-h/2}^{h/2} -y\sigma \, dy = 2b \int_{0}^{h/2} -y\sigma \, dy \tag{3.29} \\ &= 2b \left( \int_{0}^{c} -y \left(-\frac{y}{c}\sigma_0\right) dy + \int_{c}^{h/2} y\sigma_0 \, dy \right) \nonumber \end{align} $$
<figure><img src="image_53.png" /></figure>

avec :

$$ \int_{c}^{h/2} y\sigma_{0} \, dy = \sigma_{0}\left(\frac{h^{2}}{8} - \frac{c^{2}}{2}\right) \qquad (3.30) $$

En notant respectivement le moment quadratique de la zone élastique $I(c)$, et les moments statiques de la zone élastique $Z(c)$ et de la zone plastique $Z(h)$ tels que :

$$ I(c) = 2b \int_{0}^{c} y^{2} \, dy \quad ; \quad Z(c) = 2b \int_{0}^{c} y \, dy \quad ; \quad Z(h) = 2b \int_{0}^{h} y \, dy \qquad (3.31) $$

Il vient l'expression du moment fléchissant :

$$ M_{fz} = \sigma_0 \left( \frac{I(c)}{c} + Z(h) - Z(c) \right) = b\sigma_0 \left( \frac{h^2}{4} - \frac{c^2}{3} \right) \qquad (3.32) $$

La charge limite du système est donnée pour $c = 0$, toute la section est alors plastifiée. En faisant tendre $c$ vers 0 dans l'équation (3.32), le moment fléchissant dans la section médiane associée à l'état limite vaut alors :

$$ M_{fz}^{(L)} = \frac{\sigma_0 b h^2}{4} = \frac{F^{(L)} L}{4} \qquad (3.33) $$

Il vient alors la charge limite :

$$ F^{(L)} = \frac{\sigma_0 b h^2}{L} \qquad (3.34) $$

## Étude de la zone plastique

L'abscisse de la première section qui commence à plastifier est repérée par l'abscisse $a = L/2 = 2M_{fz}^{(1)}/F^{(1)}$ (ou dit autrement $M_{fz}^{(1)} = F^{(1)}L/4$, donc l'équilibre), c'est-à-dire la section médiane. À mesure que la charge croît, la zone plastique se propage aux sections voisines, de sorte que $a = 2M_{fz}/F < L/2$ avec $M_{fz}^{(1)} < M_{fz} \le M_{fz}^{(L)}$ et $F^{(1)} < F \le F^{(L)}$. La forme de cette zone est définie par la cote $c(x)$. Dans le tronçon de poutre $a \le x \le L/2$, le moment fléchissant est donné par l'équilibre :

$$ M_{fz} = \frac{Fx}{2} \qquad (3.35) $$

En combinant (3.35) et (3.32), et en y introduisant l'expression de la charge limite (3.34), il vient l'expression de $c(x)$ en fonction du rapport de charge $F/F^{(L)}$ :

$$ c(x) = h \sqrt{\frac{3}{2} \left( \frac{1}{2} - \frac{F}{F^{(L)}} \frac{x}{L} \right)} \qquad (3.36) $$

L'équation est tracée sur la figure 3.11, la zone grisée étant la zone plastifiée dans la poutre qui est limitée par une parabole définie pour $F^{(1)} < F \le F^{(L)}$. On vérifie d'autre part dans (3.36) qu'en $x = L/2$, en faisant tendre $F$ vers $F^{(L)}$, $c$ tend vers zéro.

Aussi, la zone plastique s'étend jusqu'à ce que la section médiane soit totalement plastifiée, la charge limite $F^{(L)}$ est alors atteinte : il y a écoulement libre. La propagation de la zone plastique est illustrée sur la figure 3.12.
<figure><img src="image_54.png" /><figcaption>**Figure 3.11** – Zone plastique en flexion simple</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**Figure 3.12** – Propagation de la zone plastique</figcaption></figure>

**Étude de la déformation dans la poutre** Dans la phase élastoplastique, l'étude de la déformation est conduite dans deux zones :

— la zone élastique $x \le a$, pour laquelle la loi de comportement élastique permet de relier le moment fléchissant et la courbure, autrement dit la dérivée seconde du déplacement transverse $v''(x)$ pour des poutres d'Euler-Bernoulli :

$$EIv''(x) = M_{fz}(x) \qquad (3.37)$$

Le champ de déplacement est alors obtenu par double intégration de la courbure.

— la zone élastoplastique $a \le x \le L/2$, pour laquelle les conditions de compatibilité géométriques demeurent valides. Ces conditions relient les composantes du tenseur des déformations linéarisé $\epsilon$ aux déformations et courbures de la poutre, soit considérant la déformation longitudinale :

$$\epsilon_{xx} = \frac{du_x}{dx} + \gamma_y z - \gamma_z y \qquad (3.38)$$

où $\gamma_y$ et $\gamma_z$ sont les courbures autour des directions y et z respectivement. Sous l'hypothèse d'Euler-Bernoulli, la courbure $\gamma_z$ est égale à la dérivée seconde du déplacement transverse $\gamma_z = v''(x)$. Dans notre cas, seule la flexion est à prendre en compte, elle ne s'opère d'autre part que dans le plan $(x, y)$ ce qui entraine $\epsilon_{xx} = -\gamma_z y$. La zone plastique est repérée par la cote $c(x)$ définie dans l'équation (3.36), et à cette ordonnée, le critère est atteint, soit $\sigma = -\sigma_0$. En utilisant la loi d'élasticité, toujours valable en ce point, et la condition de compatibilité, il vient :

$$\sigma(a \le x \le \frac{L}{2}, y = c) = -\sigma_0 = E\epsilon = -Ec\gamma_z \quad \text{soit} \quad c = \frac{\sigma_0}{E\gamma_z} \qquad (3.39)$$
<figure><img src="image_56.png" /></figure>

Or la courbure $\gamma_z$ peut être reliée au moment fléchissant $M_{fz}$ au travers de la cote c en introduisant l'équation (3.39) dans (3.32), il vient en fonction du moment limite (3.33) et de la courbure de début de plastification $\gamma_z^{(1)} = 2\sigma_0/(Eh)$ :

$$M_{fz} = M_{fz}^{(L)} \left( 1 - \frac{1}{3} \left( \frac{\gamma_z^{(1)}}{\gamma_z} \right)^2 \right) \qquad (3.40)$$

ce qui définit la loi de comportement dans la zone élastoplastique. On peut aussi écrire cette loi en exprimant la courbure en fonction du moment fléchissant :

$$\gamma_z = \frac{\gamma_z^{(1)}}{\sqrt{3\left(1 - \frac{M_{fz}}{M_{fz}^{(L)}}\right)}} \qquad (3.41)$$

Le champ de déplacement s'obtient alors par double intégration de (3.41). On observe dans (3.41) que la courbure $\gamma_z$ tend vers l'infini quand le moment fléchissant tend vers l'état limite $M_{fz}^{(L)}$. Aussi, tout se passe comme si la poutre était formée de deux tronçons articulés : il se forme une rotule plastique. La figure 3.13 représente l'évolution de la courbure pour les trois phases de chargement. La courbure adopte un comportement asymptotique et tend vers l'infini lorsque la charge tend vers la charge limite $F^{(L)}$.

<figure><img src="image_57.png" /><figcaption>Figure 3.13 – Évolution de la courbure au cours de la charge</figcaption></figure>

L'équation de déformée de la structure peut être calculée en intégrant les courbures (3.37) et (3.41) dans les zones respectivement élastique et élastoplastique du tronçon de poutre $0 \le x \le L/2$, en tenant compte des conditions aux limites en $x = 0, x = L/2$ et des conditions de continuité de la flèche et de la pente en $x = a$.

**Étude de la décharge** On note $\Delta F = 0 - F^{(L)}$ la variation de la charge $F$ à partir de la charge limite. Dans un premier temps, la décharge est supposée élastique. Aussi la variation du moment fléchissant dans le tronçon $0 \le x \le L/2$ s'écrit :

$$\Delta M_{fz} = \frac{\Delta F x}{2} = -\frac{F^{(L)} x}{2} \qquad (3.42)$$

de sorte qu'après la décharge supposée élastique le moment s'écrive dans ce tronçon de poutre :

$$M_{fz} = M_{fz}^{(L)} + \Delta M_{fz} = \frac{\sigma_0 b h^2}{2} \left( \frac{1}{2} - \frac{x}{L} \right) \qquad (3.43)$$
Il subsiste donc ici aussi un moment fléchissant résiduel à la décharge, et donc des contraintes résiduelles. Il reste à vérifier si la décharge est bien élastique en appliquant le critère :

$$|\sigma| - \sigma_0 \le 0, \quad \sigma = -\frac{M_{fz}y}{I} \qquad (3.44)$$

Le critère nous dit que la limite d'élasticité sera atteinte à la décharge en $x = L/2$.

**Courbe de charge** La figure 3.14 montre la courbe de charge de la poutre. Ici encore, on observe un « écrouissage de structure » durant la phase élastoplastique, provenant des zones encore non plastifiées de la poutre, avant d'atteindre la charge limite et un état d'écoulement libre. Enfin, la flèche résiduelle à décharge (dont le calcul n'est pas présenté ici) peut être

<figure><img src="image_58.png" /><figcaption>Figure 3.14 – Courbe de charge de la structure</figcaption></figure>

calculée en superposant une décharge élastique à la flèche maximale calculée par double intégration de (3.37) et (3.41) dans les zones respectivement élastique et élastoplastique du tronçon de poutre $0 \le x \le L/2$.

## 3.5 État résiduel d'un système élastoplastique

### 3.5.1 Incompatibilité des déformations plastiques et contraintes résiduelles

D'une façon générale, les déformations plastiques apparaissant dans une structure ne vérifient pas les équations de compatibilité, qui s'écrivent :

$$\pi_{ijk} \pi_{lmn} \frac{\partial^2 \varepsilon_{jm}}{\partial x_k \partial x_n} = 0 \qquad (3.45)$$

où $\pi_{ijk}$ représente le tenseur des permutations (ou de Lévi-Civita), et vaut 1 pour toutes permutations paires de (1, 2, 3), -1 pour toutes permutations impaires et 0 si au moins deux indices sont égaux. Ceci entraîne qu'il n'existe pas de champ de déplacement associé à ce champ de déformation plastique, soit formellement :

$$\nabla \mathbf{u}^P \text{ tel que } \boldsymbol{\varepsilon}^P = \frac{1}{2}(\nabla \mathbf{u}^P + \nabla^\mathrm{T} \mathbf{u}^P) \qquad (3.46)$$

On va montrer ici que l'incompatibilité des déformations plastiques entraîne l'existence de contraintes résiduelles qui subsistent après la décharge (supposée élastique) — donc un champ d'autocontrainte —, ou plutôt que la compatibilité de ces déformations impliquerait des contraintes résiduelles nulles.
<figure><img src="image_59.png" /></figure>

Considérons pour cela une structure occupant un domaine $\Omega$ constitué d'un matériau élastoplastique, soumise à des forces surfaciques données $\mathbf{F}^d$ sur tout son bord $\partial\Omega$. Le problème résolu en élasticité pure donne la solution $(\boldsymbol{\sigma}^{\text{el}}, \mathbf{u}^{\text{el}})$ :

$$
\begin{align}
\operatorname{div} \boldsymbol{\sigma}^{\text{el}} &= \mathbf{0} & \quad \forall \mathbf{x} \in \Omega \nonumber \\
\boldsymbol{\sigma}^{\text{el}} &= \mathbf{K} : \boldsymbol{\varepsilon}^{\text{el}} & \quad \forall \mathbf{x} \in \Omega \nonumber \\
\boldsymbol{\varepsilon}^{\text{el}} &= (\nabla \mathbf{u}^{\text{el}} + \nabla^{\text{T}} \mathbf{u}^{\text{el}})/2 & \quad \forall \mathbf{x} \in \Omega \tag{3.47} \\
\boldsymbol{\sigma}^{\text{el}} \cdot \mathbf{n} &= \mathbf{F}^d & \quad \forall \mathbf{x} \in \partial\Omega \nonumber
\end{align}
$$

La déformation obtenue en élasticité pure $\boldsymbol{\varepsilon}^{\text{el}}$ n'a rien à voir avec la partie élastique de la déformation $\boldsymbol{\varepsilon}^{\text{e}}$ dans le problème élastoplastique. Il en va de même pour le champ de contrainte $(\boldsymbol{\sigma}^{\text{el}} \neq \boldsymbol{\sigma})$. On définit alors les contraintes résiduelles telles que :

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\text{el}} + \boldsymbol{\sigma}^{\text{res}} \qquad (3.48) $$

La contrainte $\boldsymbol{\sigma}$ est la contrainte générée dans la structure le sous chargement $\mathbf{F}^d$ dans le problème élastoplastique. Une fois la charge effectuée, on effectue une décharge de sorte que (i) le chargement varie de $\mathbf{F}^d$ à $\mathbf{0}$ et (ii) on ajoute à $\boldsymbol{\sigma}$ les contraintes calculées par le problème purement élastique pour le chargement $-\mathbf{F}^d$, i.e. $-\boldsymbol{\sigma}^{\text{el}}$. Il vient l'état de contrainte final, résiduel $\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{el}} = \boldsymbol{\sigma}^{\text{res}}$.

Si le champ de contraintes résiduelles est nul $(\boldsymbol{\sigma}^{\text{res}} = \mathbf{0})$, il s'ensuit que $\boldsymbol{\sigma} = \boldsymbol{\sigma}^{\text{el}}$. Or ces contraintes étant reliées par la loi d'élasticité respectivement à la partie élastique des déformations $\boldsymbol{\varepsilon}^{\text{e}}$ (pour le problème élastoplastique) et à la déformation du problème purement élastique $\boldsymbol{\varepsilon}^{\text{el}}$, il s'ensuit que $\boldsymbol{\varepsilon}^{\text{e}} = \boldsymbol{\varepsilon}^{\text{el}}$. La déformation du problème purement élastique $\boldsymbol{\varepsilon}^{\text{el}}$ étant la déformation totale de ce problème, elle est compatible (elle vérifie les équations (3.45)). Il s'ensuit alors que la partie élastique des déformations $\boldsymbol{\varepsilon}^{\text{e}}$ du problème élastoplastique est compatible. Or comme la déformation totale $\boldsymbol{\varepsilon}$ est toujours compatible, par l'additivité des déformations, il vient que $\boldsymbol{\varepsilon}^{\text{p}}$ est compatible si les contraintes résiduelles sont nulles.

$$ \boldsymbol{\varepsilon}^{\text{p}} \text{ compatible} \rightarrow \boldsymbol{\sigma}^{\text{res}} = \mathbf{0} \qquad (3.49) $$

Dans ce cas, les équations vérifiées par $(\boldsymbol{\sigma}, \boldsymbol{\varepsilon}^{\text{e}}, \mathbf{u}^{\text{e}})$ sont identiques à celles vérifiées par $(\boldsymbol{\sigma}^{\text{el}}, \boldsymbol{\varepsilon}^{\text{el}}, \mathbf{u}^{\text{el}})$ (3.47), soit :

$$
\begin{align}
\operatorname{div} \boldsymbol{\sigma} &= \mathbf{0} & \quad \forall \mathbf{x} \in \Omega \nonumber \\
\boldsymbol{\sigma} &= \mathbf{K} : \boldsymbol{\varepsilon}^{\text{e}} & \quad \forall \mathbf{x} \in \Omega \nonumber \\
\boldsymbol{\varepsilon}^{\text{e}} &= (\nabla \mathbf{u}^{\text{e}} + \nabla^{\text{T}} \mathbf{u}^{\text{e}})/2 & \quad \forall \mathbf{x} \in \Omega \tag{3.50} \\
\boldsymbol{\sigma} \cdot \mathbf{n} &= \mathbf{F}^d & \quad \forall \mathbf{x} \in \partial\Omega \nonumber
\end{align}
$$

La solution du problème d'élasticité est unique, à un mouvement de corps rigide près :

$$
\left.
\begin{array}{l}
\boldsymbol{\sigma} = \boldsymbol{\sigma}^{\text{el}} \\
\boldsymbol{\varepsilon}^{\text{e}} = \boldsymbol{\varepsilon}^{\text{el}} \\
\mathbf{u}^{\text{e}} = \mathbf{u}^{\text{el}} + \text{mvt de corps rigide}
\end{array}
\right| \rightarrow \boldsymbol{\sigma}^{\text{res}} = \mathbf{0} \qquad (3.51)
$$

En général, $\boldsymbol{\varepsilon}^{\text{p}}$ est incompatible, ce qui génère des contraintes résiduelles non nulles $\boldsymbol{\sigma}^{\text{res}} \neq \mathbf{0}$. Il existe une exception : le cas où la déformation est homogène (traction simple, traction-torsion sur tube mince). La déformation $\boldsymbol{\varepsilon}^{\text{p}}$ est alors compatible et les contraintes résiduelles sont nulles $\boldsymbol{\sigma}^{\text{res}} = \mathbf{0}$.

Les différentes « déformations » constituant la partition additive de la déformation totale ne sont pas en général de véritables déformations puisqu'elles ne sont pas associées à des champs de
déplacement cinématiquement admissibles. Aussi l'incompatibilité des déformations plastiques est en tous points analogue à l'incompatibilité des déformations thermiques en thermoélasticité, des déformations hygrométriques en hygroélasticité (gonflement par l'eau), des déformations de transformation lorsqu'il y a des transformations de phase, etc. Cette incompatibilité est à l'origine de contraintes résiduelles dans les structures. Dans le cas de l'élastoplasticité, il faut ajouter aux déformations plastiques la partie élastique de la déformation, elle-même incompatible, qui rend compatible la somme des deux termes qui ne le sont pas, à la fois géométriquement et cinématiquement. Cependant, ce complément de déformation nécessaire à la compatibilité de la déformation totale génère des contraintes $\sigma = K : \varepsilon^e$. Ce champ est solution du problème homogène d'élasticité linéaire où $\varepsilon^p$ joue le rôle d'un champ de déformation initiale. Ce champ est alors identifié au champ de contraintes résiduelles $\sigma^{\text{res}}$; ces contraintes sont associées à l'élasticité du matériau. Si on néglige l'élasticité, les déformations plastiques deviennent compatibles et les contraintes résiduelles sont indéterminées.

La détermination de l'état résiduel d'une structure est d'une grande importance notamment en simulation numérique des procédés de fabrication, celui-ci détermine l'état initial de la structure lorsqu'elle sera soumise à un chargement en service. D'une part les contraintes résiduelles peuvent jouer un rôle prépondérant dans la tenue mécanique de la structure, et conduire à sa fragilisation voire à sa rupture prématurée. D'autre part, les distorsions résiduelles génèrent quant à elles des défauts géométriques qu'il convient de maîtriser pour l'assemblage de grandes structures.

L'émergence de contraintes résiduelles dans la structure est parfois recherchée. En effet, les contraintes de compression à la surface d'une structure ont tendance à limiter l'amorçage de fissures et à freiner leur propagation. C'est notamment le cas des pièces mécaniques subissant d'importants efforts de contact (éléments roulants, billes et bagues de roulements...). Ces contraintes peuvent être générées, entre autres, par grenaillage, grenaillage par choc laser ou par des traitements thermiques superficiels.

### 3.5.2 Paramètres généralisés de chargement

Dans beaucoup de problèmes, le chargement ne dépend que d'un nombre fini de paramètres, dits *paramètres généralisés de chargement* $Q = (Q_1, \dots, Q_n)$. À ces paramètres sont associés des paramètres cinématiques généralisés $q = (q_1, \dots, q_n)$, qui sont tels que la puissance des efforts extérieurs s'écrit :

$$P_{\text{ext}} = Q \cdot \dot{q} \qquad (3.52)$$

Le paramètre de chargement pourra être une force ou un couple, le paramètre cinématique associé sera alors respectivement un déplacement ou un angle.

### 3.5.3 Comportement irréversible du système

Avec ces définitions de paramètres généralisés duax en main, il est possible de décrire le comportement d'un système (donc à l'échelle de la structure). Si on considère par exemple que le système a un comportement *purement* élastique, ces paramètres peuvent être reliés tel que :

$$q(t) = \Lambda \cdot Q(t) \qquad (3.53)$$

où $\Lambda$ désigne la matrice symétrique des complaisances élastiques du système.

On envisage à présent un comportement élastoplastique du système, de sorte que $Q(t)$ soit sur la frontière actuelle du domaine d'élasticité et $\dot{Q}(t)$ soit dirigé vers l'extérieur. Le champ de
<figure><img src="image_60.png" /></figure>

taux de déformation pour une telle solution s'écrit en combinant la loi d'élasticité à la partition du taux de déformation :

$$ \dot{\boldsymbol{\epsilon}} = \boldsymbol{S} : \dot{\boldsymbol{\sigma}} + \dot{\boldsymbol{\epsilon}}^{\boldsymbol{P}} \qquad (3.54) $$

Du point de vue du comportement incrémental du système, on peut décomposer de façon analogue $\dot{\boldsymbol{q}}(t)$ en sa partie élastique donnée en dérivant (3.53) et son complément, une partie résiduelle notée $\dot{\boldsymbol{q}}^{\text{res}}(t)$ :

$$ \dot{\boldsymbol{q}}(t) = \boldsymbol{\Lambda} \cdot \dot{\boldsymbol{Q}}(t) + \dot{\boldsymbol{q}}^{\text{res}}(t) \qquad (3.55) $$

Au niveau local, on mesure le taux de déformation résiduelle $\dot{\boldsymbol{\epsilon}}^{\text{res}}$, géométriquement compatible et cinématiquement admissible avec $\dot{\boldsymbol{q}}^{\text{res}}$, et le taux de contrainte résiduelle $\dot{\boldsymbol{\sigma}}^{\text{res}}$ statiquement admissible avec $\dot{\boldsymbol{Q}}(t) - \dot{\boldsymbol{Q}}(t) = \mathbf{0}$, c'est-à-dire qu'il est autoéquilibré. Avec ces notations, on voit que le champ de taux de contrainte résiduelle $\dot{\boldsymbol{\sigma}}^{\text{res}}$ associé à $\dot{\boldsymbol{Q}}(t)$ est donné en dérivant (3.48) :

$$ \dot{\boldsymbol{\sigma}}^{\text{res}} = \dot{\boldsymbol{\sigma}} - \dot{\boldsymbol{\sigma}}^{\text{el}} \qquad (3.56) $$

où on rappelle que $\dot{\boldsymbol{\sigma}}^{\text{el}}$ dénote le taux de contrainte de la solution *purement élastique* du problème. Le champ de taux de déformation résiduelle est alors donné par :

$$ \dot{\boldsymbol{\epsilon}}^{\text{res}} = \dot{\boldsymbol{\epsilon}} - \dot{\boldsymbol{\epsilon}}^{\text{el}} = \dot{\boldsymbol{\epsilon}} - \boldsymbol{S} : \dot{\boldsymbol{\sigma}}^{\text{el}} \qquad (3.57) $$

En substituant (3.54) et (3.56) dans cette dernière équation, on obtient l'équation fondamentale :

$$ \dot{\boldsymbol{\epsilon}}^{\text{res}} = \boldsymbol{S} : \dot{\boldsymbol{\sigma}}^{\text{res}} + \dot{\boldsymbol{\epsilon}}^{\boldsymbol{P}} \qquad (3.58) $$

soit en intégrant par rapport au temps :

$$ \boldsymbol{\epsilon}^{\text{res}} = \boldsymbol{S} : (\boldsymbol{\sigma}^{\text{res}} - \boldsymbol{\sigma}_{0}) + \boldsymbol{\epsilon}^{\boldsymbol{P}} \qquad (3.59) $$

où $\boldsymbol{\sigma}^{\text{res}}$ et $\boldsymbol{\sigma}_{0}$ sont des champs autoéquilibrés. On observe alors que les déformations élastiques générées par les contraintes résiduelles, dûes à l'incompatibilité des déformations plastiques, en s'ajoutant aux déformations plastiques produisent des déformations résiduelles qui, elles, sont compatibles. Aussi ces déformations résiduelles dérivent d'un champ de déplacement

<figure><img src="image_61.png" /><figcaption>**Figure 3.15** – Déformations résiduelles : $\boldsymbol{\sigma}^{\text{res}}$ champ statiquement admissible avec $\boldsymbol{Q}(t) = \mathbf{0}$ et $\boldsymbol{\sigma}^{\text{res}}$ champ autoéquilibré</figcaption></figure>

résiduel $\boldsymbol{u}^{\text{res}}$ qu'on nomme généralement les distorsions résiduelles. Il est important d'observer que les déformations résiduelles et les contraintes résiduelles ne sont pas liées par la loi d'élasticité mais par la relation (3.59).
### 3.5.4 Écrouissage de structure en variables généralisées

On a observé à la section 3.4, en traçant les courbes de structure du treillis trois barres et de la poutre en flexion simple, un écrouissage de structure bien que considérant un matériau élastique parfaitement plastique. Un résultat général montre que les structures élastoplastiques constituées d'un matériau élastique parfaitement plastique ou à écrouissage positif soumises à un chargement, qu'on peut représenté en termes de variables généralisées (Q, q), vérifient l'inégalité suivante si le postulat de Drucker (1.49) est vérifié :

$$ \dot{\mathbf{Q}} \cdot \dot{\mathbf{q}} \geq 0 \qquad (3.60) $$

Ce résultat se montre facilement en appliquant le principe des puissances virtuelles appliqué au champ réel :

$$ \int_{\Omega} \sigma : \dot{\mathbf{e}} \, d\Omega = \mathbf{Q} \cdot \dot{\mathbf{q}} \qquad (3.61) $$

en différenciant par rapport au temps et en utilisant (3.54), il vient :

$$ \dot{\mathbf{Q}} \cdot \dot{\mathbf{q}} = \int_{\Omega} \dot{\sigma} : \dot{\mathbf{e}} \, d\Omega = \int_{\Omega} \dot{\sigma} : \mathbf{S} : \dot{\sigma} \, d\Omega + \int_{\Omega} \dot{\sigma} : \mathbf{\epsilon}^{\mathrm{P}} \, d\Omega \qquad (3.62) $$

Les deux premiers termes de cette égalité sont séparément positifs, le premier en raison de la définie-positivité de S, le second en vertu du postulat de Drucker (1.49) établi au premier chapitre.

La relation (3.62) permet d'expliquer l'origine de l'écrouissage visible sur la courbe de structure (q, Q) :

1. la première contribution est celle de l'écrouissage propre du matériau. Cette contribution intervient dans le second terme de (3.62). Comme expliqué à la section 1.5.4, $\dot{\sigma} : \dot{\mathbf{e}}^{\mathrm{P}}$ est strictement positif lorsque l'écrouissage du matériau est aussi strictement positif, et nul pour un matériau élastique parfaitement plastique;

2. la seconde contribution a trait à ce qu'on appelle « l'écrouissage de structure » représenté par le terme $\dot{\sigma} : \mathbf{S} : \dot{\sigma}$ dans l'équation (3.62). Ce terme est strictement positif tant qu'il y a évolution des contraintes, i.e. $\dot{\sigma} \neq 0$. Cette contribution est dûe :

- à la présence d'une zone élastique au sein de la structure, permettant aux contraintes d'y croître. Cette croissance s'accompagne aussi d'une extension de la zone plastique;

- elle peut aussi être due à la rotation des contraintes sur le seuil de plasticité, ceci peut apparaître lorsque les équations d'équilibre ne déterminent pas complètement l'état de contrainte.
Simulation numérique des
structures élastoplastiques

4.1 Simulation numérique en calcul des structures

De manière générale, l'amélioration des capacités de calcul et des outils numériques permet d'aborder des problèmes de plus en plus compliqués, et de réaliser des simulations de difficulté croissante. Les quelques solutions analytiques connues de problèmes de l'ingénieur se limitent à des cas « académiques ». Les *problèmes industriels*, quant à eux, impliquent des géométries complexes, des comportements non linéaires avec effets d'histoire, des structures multi-matériaux et/ou multi-échelles, de la fissuration, etc. La solution repose alors nécessairement sur la simulation numérique.

La simulation numérique est inégalée pour la bonne compréhension des mécanismes physiques intervenant dans une analyse, pour la conception d'un produit ou la mise au point d'un procédé ; elle facilite l'accès à des quantités mécaniques non mesurables par l'expérience, comme, par exemple, les contraintes résiduelles. Plus généralement, la simulation numérique remplace a bien moindre coût le prototypage.

L'édification d'un modèle pose toujours la question de sa mission. La bonne recette combine ingrédients variés et techniques numériques sophistiquées : son degré de finesse est subordonné à une claire définition des objectifs. La complexité et le raffinement du modèle sont généralement mis en regard d'un critère informations/coût de calcul, ce dernier limitant son encombrement numérique.

On s'intéresse dans ce chapitre à la simulation numérique des structures élastoplastiques. Le comportement élastoplastique étant en partie irréversible, les effets d'histoire pilotent l'état de l'élément de matière et le problème élastoplastique se formule comme un problème d'évolution non linéaire pour lequel une stratégie de résolution incrémentale est nécessaire. La section 4.2 présente brièvement la résolution de tels problèmes dans le cadre de la méthode des éléments finis. La section 4.3 se concentre sur les algorithmes d'intégration des lois de comportement élastoplastiques.

4.2 Méthode des éléments finis en mécanique non-linéaire

4.2.1 Formulation faible du problème élastoplastique

Le problème d'évolution élastoplastique (P) peut se réécrire sous forme faible de la façon suivante : Trouver le couple déplacement-contraintes (u(x, t), σ(x, t)) défini en tout point et tout instant (x, t) ∈ Ω × [0, T], appartenant à Uad[0,T] × Sad[0,T] et solution :
— des équations de liaisons et conditions initiales :

$$
\begin{align}
\mathbf{u}(\mathbf{x}, t) &\in \mathcal{U}_{\text{ad}}^{[0,T]} \nonumber \\
\mathbf{u}(\mathbf{x}, t) &= \mathbf{u}^{\text{d}}(\mathbf{x}, t) \quad \forall \mathbf{x} \in \partial\Omega_{\text{u}}, \forall t \nonumber \\
\mathbf{u}(\mathbf{x}, t=0) &= \mathbf{u}_{0}(\mathbf{x}) \quad \forall \mathbf{x} \in \Omega \tag{4.1}
\end{align}
$$

— du Principe des Puissances Virtuelles (PPV) : trouver $\sigma \in S_{\text{ad}}^{[0,T]}$ tel que

$$
\forall \mathbf{u}^* \in \mathcal{U}_{\text{ad}}^0, \quad \int_{\Omega} \sigma : \varepsilon(\mathbf{u}^*) \, d\Omega = \int_{\Omega} \mathbf{f}_d \cdot \mathbf{u}^* \, d\Omega + \int_{\partial\Omega_F} \mathbf{F}_d \cdot \mathbf{u}^* \, dS \quad (4.2)
$$

— de la relation de comportement

$$
\sigma(\mathbf{x}, t) = \mathcal{F}(\varepsilon(\dot{\mathbf{u}}(\mathbf{x}, t)), \tau \le t) \quad (4.3)
$$

où la fonctionnelle $\mathcal{F}$ symbolise la relation incrémentale du comportement élastoplas-
tique.

Les espaces fonctionnels sont définis comme suit :

— espace des déplacements admissibles homogènes

$$
\mathcal{U}_{\text{ad}}^{0} = \{ \mathbf{u}^{*} \text{ régulier sur } \Omega \mid \mathbf{u}^{*} = \mathbf{0}, \forall \mathbf{x} \in \partial\Omega_{u}, \forall t \in [0, T] \}
$$

— espace des contraintes admissibles

$$
S_{\text{ad}} = \{ \boldsymbol{\sigma}^* \text{ régulier sur } \Omega, \sigma_{ij}^* = \sigma_{ji}^*, \sigma_{ij,j}^* + f_i = 0 \mid \boldsymbol{\sigma}^* \cdot \mathbf{n} = \mathbf{F}_d, \forall \mathbf{x} \in \partial\Omega_F, \forall t \in [0, T] \}
$$

Les équations de l'élastoplasticité étant incrémentales, la solution à un instant *t* ne peut être obtenue que par une méthode pas à pas. Considérons l'intervalle de temps [0, *T*], découpé en *M* pas de temps (supposés égaux pour simplifier) *Δt* = *T*/*M*. À l'instant *t*<sub>n</sub>, toutes les grandeurs mécaniques du problème sont connues :

$$
S_n \stackrel{\text{def}}{=} \{ \mathbf{u}_n, \boldsymbol{\epsilon}_n, \boldsymbol{\epsilon}_n^p, \boldsymbol{\sigma}_n \dots \} \tag{4.4}
$$

L'approche incrémentale consiste à déterminer $S_0, S_1, S_2, \dots, S_M$ de proche en proche. Pour un incrément, le problème s'écrit :

Connaissant $S_n$ et $(\mathbf{f}_{n+1}^d, \mathbf{u}_{n+1}^d, \mathbf{F}_{n+1}^d)$
Trouver $S_{n+1}$

La discrétisation temporelle permet de construire une succession de problèmes non linéaires et indépendants du temps. Les conditions d'admissibilité statique et cinématique sont uniquement vérifiées aux extrémités des incréments temporels. Par contre, la loi de comportement est vérifiée sur tout l'intervalle en supposant généralement une évolution linéaire des déplacements sur le pas de temps.

4.2.2 Discrétisation éléments finis et équations discrètes

On omet par la suite (sauf indication explicite) les indices temporels, le problème s'écrivant
à l'instant $t_{n+1}$, qu'on écrit aussi l'instant $t + \Delta t$. La discrétisation par éléments finis de la
forme faible (4.2) à l'instant $t + \Delta t$ conduit à l'écriture des équations discrètes suivantes :

$$
\mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}(\mathbf{u}) = \mathbf{0} \tag{4.5}
$$
<figure><img src="image_62.png" /></figure>

où **f**<sup>ext</sup> et **f**<sup>int</sup> dénotent respectivement les forces externes et internes, définies à partir des quantités élémentaires telles que :

$$ \mathbf{f}^{\text{ext}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{N}^{\top} \mathbf{f}_d \, d\Omega + \int_{\partial\Omega_e \cap \partial\Omega_F} \mathbf{N}^{\top} \mathbf{F}_d \, dS \right) \quad ; \quad \mathbf{f}^{\text{int}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{B}^{\top} \boldsymbol{\sigma} \, d\Omega \right) \quad (4.6) $$

où la notation $\sum_{e=1}^{N_e}$ dénote l'opération d'assemblage des quantités élémentaires. Les matrices **N** et **B** contiennent respectivement les fonctions de forme et les gradients des fonctions de forme de l'élément fini considéré. La résolution du système non linéaire (C.1) est effectuée sur le champ de déplacement **u**, inconnues définies aux nœuds du maillage, intervenant par l'intermédiaire de la loi de comportement élastoplastique dans les forces internes (4.6). Ce type de problème est généralement résolu par une méthode itérative, consistant à construire une suite de vecteur **u** convergente vers la solution du problème, chaque élément de cette suite étant obtenu par la résolution d'un problème linéaire.

### 4.2.3 Résolution du problème stationnaire

Le système d'équation (C.1) étant non linéaire, la résolution se fait par une méthode itérative. On écrit tout d'abord le résidu du système d'équation (C.1) :

$$ \mathbf{R}(\mathbf{u}) = \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}(\mathbf{u}) \quad (4.7) $$

En se donnant **u**<sup>(k)</sup> une solution approchée à l'itération k, la détermination d'une meilleure approximation **u**<sup>(k+1)</sup> s'obtient en effectuant une linéarisation du résidu à l'itération k du calcul avec la méthode de Newton-Raphson :

$$ \mathbf{R}^{(k+1)} \approx \mathbf{R}^{(k)} + \left. \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right|^{(k)} \delta \mathbf{u}^{(k)} = \mathbf{0} \Leftrightarrow \mathbf{R}(\mathbf{u}^{(k+1)}) \approx \mathbf{R}(\mathbf{u}^{(k)}) + \left. \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right| (\mathbf{u}^{(k)}) \delta \mathbf{u}^{(k)} = \mathbf{0} \quad (4.8) $$

où $\delta \mathbf{u}^{(k)}$ est un incrément de la solution entre deux itérations du calcul, défini tel que $\delta \mathbf{u}^{(k)} = \mathbf{u}^{(k+1)} - \mathbf{u}^{(k)}$. On résout alors à chaque itération du calcul le système linéaire suivant :

$$ \mathbf{K}^{(k)} \delta \mathbf{u}^{(k)} = \mathbf{R}^{(k)} \quad \text{avec} \quad \mathbf{K}^{(k)} = -\left. \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right|^{(k)} \quad \text{et} \quad \mathbf{R}^{(k)} = \mathbf{f}^{\text{ext}} - (\mathbf{f}^{\text{int}})^{(k)} \quad (4.9) $$

où $\mathbf{K}^{(k)}$ est la matrice tangente calculée à l'itération k du calcul, cohérente avec l'algorithme de Newton-Raphson. L'initialisation de l'algorithme (itération k = 0) est généralement effectuée en supposant un comportement élastique linéaire du matériau en tous points de la structure, la matrice élastique s'écrivant :

$$ \mathbf{K}_{\text{elas}} = \int_{\Omega} \mathbf{B}^{\top} \mathbf{D} \mathbf{B} \, d\Omega \quad (4.10) $$

où **D** est l'opérateur d'élasticité. Le couple solution ($\mathbf{u}^{(0)}$, $\boldsymbol{\sigma}^{(0)}$) obtenu est cinématiquement et statiquement admissible au sens des éléments finis, en revanche la relation de comportement n'est, elle, pas forcément vérifiée. Pour tester la qualité de la solution, il faut évaluer la valeur du résidu $\mathbf{R}^{(0)}$, et comparer sa norme à une valeur donnée par l'utilisateur :

1. la première étape est la détermination de l'état de contrainte vérifiant la loi de comportement associée à l'accroissement de déplacement calculé; cette étape correspond à l'intégration de la loi de comportement :

$$ \hat{\boldsymbol{\sigma}}^{(0)} = \boldsymbol{\sigma}_n + \mathcal{F}\left(\mathbf{B}\frac{(\mathbf{u}^{(0)} - \mathbf{u}_n)}{t_{n+1} - t_n}; \tau \in [t_n, t_{n+1}]\right) \quad (4.11) $$
2. la seconde étape est l'évaluation du résidu (4.7) tenant compte de la contrainte (4.11). La valeur du résidu caractérise la vérification des équations d'équilibre, et mesure en quelque sorte le « déséquilibre ».

Si la norme du résidu $||\mathbf{R}^{(0)}||$ est supérieure à une valeur fixée par l'utilisateur, il faut déterminer une nouvelle approximation $\mathbf{u}^{(1)}$ de la solution, en assemblant la matrice tangente et en résolvant le système linéaire (C.4). Compte tenu de l'hypothèse des petites perturbations, la matrice tangente s'écrit :

$$ \mathbf{K}^{(k)} = \int_{\Omega} \mathbf{B}^{\top} \frac{\partial \boldsymbol{\sigma}}{\partial \mathbf{u}} \bigg|^{(k)} d\Omega \quad (4.12) $$

où seule la partie « comportement » est dérivée. Au voisinage du point intermédiaire $(\boldsymbol{\varepsilon}^{(k)}, \hat{\boldsymbol{\sigma}}^{(k)})$, on peut extraire de la loi de comportement le module tangent discret $\mathbf{H}^{(k)}$, ou opérateur tangent, de sorte que :

$$ \delta\boldsymbol{\sigma}^{(k)} = \mathbf{H}^{(k)}\delta\boldsymbol{\varepsilon}^{(k)} \quad \text{avec} \quad \delta\boldsymbol{\sigma}^{(k)} = \boldsymbol{\sigma}^{(k+1)} - \hat{\boldsymbol{\sigma}}^{(k)} \quad \text{et} \quad \delta\boldsymbol{\varepsilon}^{(k)} = \mathbf{B}(\mathbf{u}^{(k+1)} - \mathbf{u}^{(k)}) \quad (4.13) $$

La matrice tangente s'écrit alors en introduisant l'opérateur tangent discret, cohérent :

$$ \frac{\partial \boldsymbol{\sigma}}{\partial \mathbf{u}} \bigg|^{(k)} = \mathbf{H}^{(k)}\mathbf{B} \rightarrow \mathbf{K}^{(k)} = \int_{\Omega} \mathbf{B}^{\top} \mathbf{H}^{(k)} \mathbf{B} \, d\Omega \quad (4.14) $$

La solution intermédiaire du champ de contraintes statiquement admissible au sens des éléments finis est donnée par :

$$ \boldsymbol{\sigma}^{(k+1)} = \hat{\boldsymbol{\sigma}}^{(k)} + \mathbf{H}^{(k)} \mathbf{B} \delta \mathbf{u}^{(k)} \quad (4.15) $$

Cependant, le couple $(\mathbf{u}^{(k+1)}, \boldsymbol{\sigma}^{(k+1)})$ ne vérifie pas la loi de comportement. On préfère en

<figure><img src="image_63.png" /><figcaption>Figure 4.1 – Résolution du problème élastoplastique par la méthode des éléments finis</figcaption></figure>

général retenir une approximation du champ de contrainte vérifiant la loi de comportement mais n'étant pas statiquement admissible au sens des éléments finis :

$$ \hat{\boldsymbol{\sigma}}^{(k+1)} = \boldsymbol{\sigma}_n + \mathcal{F} \left( \mathbf{B} \frac{\mathbf{u}^{(k+1)} - \mathbf{u}_n}{t_{n+1} - t_n} ; \tau \in [t_n, t_{n+1}] \right) \quad (4.16) $$

En résumé, chaque itération de calcul est composée de deux étapes :
<figure><img src="image_64.png" /></figure>

**Étape globale** Il s'agit de l'étape lors de laquelle les déplacements sont calculés aux nœuds du maillage de la structure, en résolvant le système linéaire (C.4). À l'itération 0, il est supposé que l'incrément des déformations plastiques $\Delta \boldsymbol{\varepsilon}^{(0)}$ est nul (calcul *purement* élastique). Dans les itérations suivantes, la valeur trouvée à l'itération précédente joue le rôle d'une déformation initiale. Les équations de la plasticité ne sont pas utilisées durant cette étape.

**Étape locale** Elle consiste à intégrer la loi de comportement en chaque point de Gauss de la structure. Connaissant l'incrément des déformations totales $\Delta \boldsymbol{\varepsilon}$, on calcule l'incrément des déformations plastiques $\Delta \boldsymbol{\varepsilon}^P$ et des variables internes $\Delta \boldsymbol{V}_k$, ainsi que les contraintes $\boldsymbol{\sigma}$ et forces thermodynamiques associées $\boldsymbol{A}_k$ en utilisant les équations de la plasticité.

Le schéma 4.1 résume le processus itératif de la résolution du problème élastoplastique par la méthode des éléments finis. L'étape globale permet d'avoir accès à une meilleure approximation $\mathbf{u}^{(k+1)}$ du champ solution, au travers de la matrice tangente (C.4) qui définit la direction de recherche à l'itération $k$ du calcul. L'évaluation du résidu effectuée avec l'état de contrainte associé à l'approximation $\mathbf{u}^{(k+1)}$ par la loi de comportement définit la « distance » entre le niveau des forces externes et internes. La figure 4.2 illustre ce processus itératif pour un système à un degré de liberté. En pratique, l'algorithme C.1 montre le processus itératif de résolution du

<figure><img src="image_65.png" /><figcaption>**Figure 4.2** – Convergence de l’algorithme pour un système à un degré de liberté</figcaption></figure>

problème d’évolution élastoplastique par la méthode des éléments finis.

1. **Données** : maillage éléments finis, historique de chargement, conditions aux limites, loi de comportement

2. **for** $n = 1$ **to** $M$ **pas de chargement** **do**

3. **prédiction élastique** : $\mathbf{u}^{(0)}$ solution of $K_{\text{elas}}\mathbf{u} = \mathbf{f}^{\text{ext}}$

4. **initialisation du résidu** : $\mathbf{R}^{(0)} \leftarrow \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}(\mathbf{u}^{(0)})$

5. **while** $\|\mathbf{R}^{(k)}\| > \text{TOL}$ **do**

6. **étape globale** : $K^{(k)}\delta\mathbf{u}^{(k)} = \mathbf{R}^{(k)} \rightarrow \mathbf{u}^{(k+1)}$

7. **for** $e = 1$ **to** $N_e$ **éléments (calcul du résidu)** **do**

8. **étape locale** : intégration de la loi de comportement

9. calcul des quantités élémentaires $K^e(\mathbf{u}^{(k+1)}), R^e(\mathbf{u}^{(k+1)}) = \mathbf{f}^e_{\text{ext}} - \mathbf{f}^e_{\text{int}}(\mathbf{u}^{(k+1)})$

10. assemblage des quantités élémentaires : $\mathbf{K} = \sum_e K^e, \mathbf{R} = \sum_e R^e$

11. **end for**

12. $k \leftarrow k + 1$

13. **end while**

14. **end for**

<caption>**Algorithme 4.1** – Résolution du problème d'évolution élastoplastique par éléments finis</caption>
## 4.3 Intégration des lois de comportement élastoplastiques

L'intégration de la loi de comportement effectuée aux points de Gauss de la structure, ou étape locale d'un calcul éléments finis, consiste à calculer l'état mécanique {$\boldsymbol{\sigma}_{n+1}, \boldsymbol{\varepsilon}_{n+1}^p$; ($\mathbf{A}_k$)$_{n+1}$, ($\mathbf{V}_k$)$_{n+1}$} à l'instant $t_{n+1}$, à partir d'un incrément de la déformation totale $\Delta \boldsymbol{\varepsilon}$ donné, et de l'état mécanique {$\boldsymbol{\sigma}_n, \boldsymbol{\varepsilon}_n^p$; ($\mathbf{A}_k$)$_n$, ($\mathbf{V}_k$)$_n$} à l'instant $t_n$.

Le processus défini sur un pas de temps $[t_n, t_{n+1}]$, est vu comme piloté par la déformation, l'incrément de déformation donné $\Delta \boldsymbol{\varepsilon}$ est calculé aux points de Gauss de l'élément à partir du champ de déplacement obtenu lors de l'étape globale, c'est-à-dire en résolvant (C.4).

<figure><img src="image_66.png" /><figcaption>**Figure 4.3** – Intégration de la loi de comportement élastoplastique en un point de Gauss. Les contraintes $\boldsymbol{\sigma}$ et les forces thermodynamiques $\mathbf{A}_k$ sont des variables dépendantes de la déformation plastique $\boldsymbol{\varepsilon}^p$ et des variables internes $\mathbf{V}_k$ au travers des lois d'état.</figcaption></figure>

### 4.3.1 Modèle de plasticité de von Mises avec écrouissage isotrope linéaire

Afin de décrire les algorithmes d'intégration des lois de comportement les plus connus par la suite, on va s'appuyer sur un modèle simple de plasticité qui va servir de fil rouge tout au long de cette section : le modèle de plasticité de von Mises avec un écrouissage isotrope linéaire, dont on rappelle ici brièvement les équations constitutives :

$$ \dot{\boldsymbol{\varepsilon}} = \dot{\boldsymbol{\varepsilon}}^e + \dot{\boldsymbol{\varepsilon}}^p \qquad (4.17a) $$

$$ \dot{\boldsymbol{\sigma}} = \lambda(\text{tr} \, \dot{\boldsymbol{\varepsilon}}^e)\mathbf{1} + 2\mu\dot{\boldsymbol{\varepsilon}}^e \qquad (4.17b) $$

$$ f(\boldsymbol{\sigma}, R) = \sqrt{\frac{3}{2}\boldsymbol{s} : \boldsymbol{s} - hp - \sigma_y} \le 0 \qquad (4.17c) $$

$$ p = \int_{0}^{T} \dot{p}(\tau) \, d\tau \quad ; \quad \dot{p} = \sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^{p} : \dot{\boldsymbol{\varepsilon}}^{p}} \qquad (4.17d) $$

$$ \dot{\boldsymbol{\varepsilon}}^{p} = \dot{p} \frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{3}{2} \dot{p} \frac{\boldsymbol{s}}{\sigma_{eq}} \qquad (4.17e) $$

$$ \dot{p} \ge 0 \quad ; \quad f \le 0 \quad ; \quad f \dot{p} = 0 \qquad (4.17f) $$

comprenant la partition additive du taux de déformation (4.17a), la loi d'élasticité (4.17b) avec $\lambda$ et $\mu$ les coefficients de Lamé, le critère de von Mises (4.17c) prenant en compte un écrouissage linéaire isotrope de paramètre $h$, la définition de la déformation plastique cumulée (4.17d), la loi d'écoulement plastique (4.17e) et enfin les conditions de complémentarité de Kühn-Tucker (4.17f).

### 4.3.2 Algorithmes de type Return Mapping

Les algorithmes de « Return Mapping » forment une classe de méthodes d'intégration robuste et précise, très utilisée en pratique, qui permet d'effectuer l'intégration numérique d'une loi de
<figure><img src="image_67.png" /></figure>

comportement à seuil. Les schémas d'intégration de type Return Mapping sont des algorithmes utilisant une étape de prédiction élastique, impliquant une excursion de l'état de contrainte hors du domaine d'élasticité, et d'une étape de correction plastique qui projette les contraintes sur la surface de charge actualisée. Ce type de méthodes est constitué d'un schéma d'intégration qui transforme un système d'équations de comportement en un système d'équations non linéaires algébriques et d'un schéma de résolution de ce système d'équations. Les équations de la plasticité (4.17) forment un système d'équations différentielles ordinaires non linéaire. La discrétisation, ou la forme incrémentale, d'un tel système peut être effectuée par les schémas classiques comme la méthode des trapèzes généralisés, des différences généralisées ou les méthodes de Runge-Kutta.

### 4.3.2.1 Retour radial

La méthode du retour radial pour la plasticité de von Mises est un cas particulier de la classe des méthodes de type Return Mapping et correspond à une discrétisation temporelle de type Euler implicite. Le système (4.17) discrétisé avec ce schéma s'écrit :

$$ \Delta \boldsymbol{\varepsilon} = \Delta \boldsymbol{\varepsilon}^e + \Delta \boldsymbol{\varepsilon}^p \quad (4.18a) $$

$$ \Delta \boldsymbol{\sigma} = \lambda(\mathrm{tr} \, \Delta \boldsymbol{\varepsilon}^e) \mathbf{1} + 2\mu \Delta \boldsymbol{\varepsilon}^e \quad (4.18b) $$

$$ f_{n+1} = \sigma_{eq_{n+1}} - hp_{n+1} - \sigma_y \le 0 \quad (4.18c) $$

$$ p_{n+1} = p_n + \int_{t_n}^{t_{n+1}} \dot{p}(\tau) \, d\tau \quad ; \quad \Delta p = \sqrt{\frac{2}{3} \Delta \boldsymbol{\varepsilon}^p : \Delta \boldsymbol{\varepsilon}^p} \quad (4.18d) $$

$$ \Delta \boldsymbol{\varepsilon}^p = \Delta p \left. \frac{\partial f}{\partial \boldsymbol{\sigma}} \right|_{n+1} = \frac{3}{2} \Delta p \frac{s_{n+1}}{\sigma_{eq_{n+1}}} \quad (4.18e) $$

$$ \Delta p \ge 0 \quad ; \quad f_{n+1} \Delta p = 0 \quad (4.18f) $$

où on rappelle que $\Delta \boldsymbol{\varepsilon}$ est une donnée. La solution du problème (4.18) passe d'abord par l'introduction d'un problème auxiliaire : l'état élastique « prédit ». Cet état ne correspond pas forcément à l'état réel et est obtenu en « gelant » l'écoulement plastique et l'évolution des variables internes. Autrement dit, on effectue un test, les quantités « test » seront affectées de la notation $(\bullet)^*$; le test associé à la prédiction purement élastique est défini comme suit :

$$ \sigma_{n+1}^* = \sigma_n + \lambda(\mathrm{tr} \, \Delta \boldsymbol{\varepsilon}) \mathbf{1} + 2\mu \Delta \boldsymbol{\varepsilon} \quad (4.19a) $$

$$ \boldsymbol{\varepsilon}_{n+1}^{p*} = \boldsymbol{\varepsilon}_n^p \quad (4.19b) $$

$$ (\mathbf{V}_k)_*_{n+1} = (\mathbf{V}_k)_*^n \quad (4.19c) $$

$$ f_{n+1}^* = \sigma_{eq_{n+1}}^* - hp_n - \sigma_y \quad (4.19d) $$

Cet état prédit ne dépend que des conditions initiales $\{\boldsymbol{\varepsilon}_n, \boldsymbol{\varepsilon}_n^p, (\mathbf{V}_k)_n\}$ à l'instant $t_n$ et de l'incrément de déformation $\Delta \boldsymbol{\varepsilon}$ donné.

Une fois l'état prédit calculé, cette prédiction est testée sur le critère de plasticité. Si cet état est admissible au sens du critère, c'est-à-dire $f_{n+1}^* \le 0$, les variables sont actualisées avec les valeurs prédites :

$$ \begin{aligned} \sigma_{n+1} &= \sigma_{n+1}^* \\ \boldsymbol{\varepsilon}_{n+1}^p &= \boldsymbol{\varepsilon}_{n+1}^{p*} \\ (\mathbf{V}_k)_{n+1} &= (\mathbf{V}_k)_*^n \end{aligned} \quad (4.20) $$
cet état satisfaisant les lois d'état, la loi d'écoulement plastique et les lois d'écrouissage car $\Delta p = 0$, et les conditions de Kühn-Tucker car $f_{n+1} = f_{n+1}^* \le 0$ et $\Delta p = 0$. L'état prédit est alors la solution du problème.

Considérons à présent le cas où l'état prédit n'est pas admissible au sens du critère $f_{n+1}^* > 0$. La solution prédite ne peut alors pas être une solution du problème incrémental. On effectue alors une étape de correction plastique dans laquelle l'objectif est de calculer la solution $\{\boldsymbol{\epsilon}_{n+1}^P, (V_k)_{n+1}, \boldsymbol{\sigma}_{n+1}, \Delta p\}$. On exprime tout d'abord l'état de contrainte final en fonction de l'état de contrainte prédit :

$$ 
\begin{align} 
\boldsymbol{\sigma}_{n+1} &= \boldsymbol{\sigma}_n + \lambda(\operatorname{tr}(\Delta\boldsymbol{\epsilon} - \Delta\boldsymbol{\epsilon}^P))\mathbf{1} + 2\mu(\Delta\boldsymbol{\epsilon} - \Delta\boldsymbol{\epsilon}^P) \nonumber \\
&= \boldsymbol{\sigma}_n + \lambda(\operatorname{tr}(\Delta\boldsymbol{\epsilon})\mathbf{1} + 2\mu\Delta\boldsymbol{\epsilon} - 2\mu\Delta\boldsymbol{\epsilon}^P) \tag{4.21} \\
&= \boldsymbol{\sigma}_{n+1}^* - 2\mu\Delta\boldsymbol{\epsilon}^P \nonumber
\end{align} 
$$

où on a utilisé l'incompressibilité des déformations plastiques $\operatorname{tr} \Delta\boldsymbol{\epsilon}^P = 0$ induite par le critère de von Mises. La quantité $2\mu\Delta\boldsymbol{\epsilon}^P$ intervient ici comme un correcteur plastique, c'est la correction qui doit être apportée à la prédiction élastique pour atteindre l'état de contrainte final. La contrainte élastique prédite est projetée sur la surface de charge finale avec le correcteur plastique ; d'autre part cette projection est normale à la surface de charge finale inconnue. Cette projection est illustrée sur la figure 4.4. Combinons la partie déviatorique de (4.21) et la loi

<figure><img src="image_68.png" /><figcaption>Figure 4.4 – Retour radial</figcaption></figure>

d'écoulement plastique (4.18e) :

$$ s_{n+1} = s_{n+1}^* - 2\mu\Delta\boldsymbol{\epsilon}^P = s_{n+1}^* - 3\mu\frac{\Delta p s_{n+1}}{\sigma_{eq_{n+1}}} \rightarrow \left(1 + \frac{3\mu\Delta p}{\sigma_{eq_{n+1}}}\right)s_{n+1} = s_{n+1}^* \quad (4.22) $$

L'équation (4.22) illustre une caractéristique fondamentale d'un algorithme implicite de projection sur la surface de charge de von Mises : le déviateur des contraintes à l'instant $t_{n+1}$ est positivement colinéaire au prédicteur élastique $s_{n+1}^*$. C'est l'algorithme de retour radial. Cette propriété simplifie grandement la construction de la solution et n'est valable qu'avec le critère de von Mises et une discrétisation implicite.
Prenons à présent la norme de von Mises de l'équation (4.22), il vient :

$$
\begin{gather}
\left(1 + 3\mu \frac{\Delta p}{\sigma_{\text{eq}n+1}}\right) \sigma_{\text{eq}n+1} = \sigma_{\text{eq}n+1} + 3\mu \Delta p = \sigma_{\text{eq}n+1}^* \tag{4.23} \\
\sigma_{\text{eq}n+1}^* = \sqrt{\frac{3}{2} s_{n+1}^* : s_{n+1}^*} \nonumber
\end{gather}
$$

Du critère (4.18c), supposé atteint, vient l’expression de l’incrément de déformation plastique cumulée :

$$
\Delta p = \frac{f_{n+1}^*}{h + 3\mu} \tag{4.24}
$$

La simplicité du modèle envisagé permet de déduire explicitement l’expression de l’incrément de déformation plastique cumulée $\Delta p$, en particulier le choix d’un écrouissage linéaire. Dans le cas plus général d’un écrouissage non linéaire, l’équation non linéaire suivante doit être résolue :

$$
\sigma_{\text{eq}n+1}^* - 3\mu\Delta p - R_{n+1}(p) - \sigma_y = 0 \quad (4.25)
$$

où $R_{n+1}(p)$ représente l'évolution de l'écrouissage isotrope évaluée à l'instant $t_{n+1}$. Tout
le problème de l'intégration numérique de la loi de comportement élastoplastique revient à
résoudre l'équation (4.25) sur $\Delta p$.

Une fois $\Delta p$ connu, l'incrément de déformation plastique est calculé avec la loi d'écoulement plastique (4.18e) combinée avec la colinéarité positive de $\mathbf{s}_{n+1}^*$ et $\mathbf{s}_{n+1}$ :

$$
\Delta \boldsymbol{\epsilon}^P = \frac{3}{2} \frac{\Delta p}{\sigma_{eq,n+1}^*} \mathbf{s}_{n+1}^* \quad (4.26)
$$

L'état de contrainte final est ensuite calculé par (4.21).

**Calcul de l'opérateur tangent cohérent** La convergence de l’algorithme itératif est grande-
ment améliorée si on utilise l’opérateur tangent élastoplastique plutôt que la matrice élastique.
Le calcul de cet opérateur tangent est détaillé pour le cas simple envisagé.

Pour une loi de comportement donnée, il est possible de déterminer l'opérateur tangent « continu » $\mathbf{H}_c$ vérifiant :

$$
\dot{\boldsymbol{\sigma}} = \boldsymbol{H}_c : \dot{\boldsymbol{\epsilon}} \tag{4.27}
$$

lorsque le couple $(\boldsymbol{\epsilon}, \boldsymbol{\sigma})$ satisfait la loi de comportement. Comme la loi de comportement est non linéaire, son intégration sur l'incrément est faite numériquement. La solution obtenue est donc approchée et dépend du schéma d'intégration utilisé. À une itération $k$ du calcul, l'opérateur tangent discret $\boldsymbol{s}^{(k)}$ s'écrit :

$$
\delta \boldsymbol{\sigma}^{(k)} = \mathbf{H}^{(k)} \delta \boldsymbol{\varepsilon}^{(k)} \quad (4.28)
$$

Les opérateurs tangents discret et continu sont identiques dans la limite où $\Delta t \to 0$. En pratique, le pas de temps $\Delta t$ n'est pas infiniment petit et la différence entre les deux opérateurs peut être importante. Il a été montré que l'opérateur discret permet de conserver la convergence quadratique de la méthode de Newton. Si l'incrément « plastique » est important, l'utilisation de l'opérateur continu à la place de l'opérateur discret détériore considérablement les performances de convergence de l'algorithme de résolution globale. L'incrément de déformation $\Delta \boldsymbol{\epsilon}$ sert
d'entrée pour l'intégration de la loi de comportement et on écrit l'opérateur tangent discret, en
omettant l'indice *k* du calcul, comme suit :

$$
\mathbf{H} = \frac{\partial \Delta \boldsymbol{\sigma}}{\partial \Delta \boldsymbol{\epsilon}} \qquad (4.29)
$$

Partant de (4.21) en appliquant la définition de l'opérateur discret (4.29), il vient :

$$
\mathbf{H} = \mathbf{D} - 2\mu \frac{\partial \Delta \boldsymbol{\epsilon}^p}{\partial \Delta \boldsymbol{\epsilon}} \qquad (4.30)
$$

où **D** est l'opérateur d'élasticité. La loi d'écoulement plastique (4.26) injectée dans (4.30)
implique :

$$
\mathbf{H} = \mathbf{D} - 3\mu \frac{\partial}{\partial \Delta \boldsymbol{\epsilon}} \left( \Delta p \frac{\mathbf{s}_{n+1}^*}{\sigma_{\mathrm{eq}\,n+1}^*} \right) = \mathbf{D} - \mathbf{E} \quad (4.31)
$$

où **E** apparaît comme « correction plastique », tenseur d'ordre 4 qu'il faut à présent calculer. Ce calcul fait intervenir le calcul préliminaire de trois quantités :

1. Calcul de $\frac{\partial s_{n+1}^*}{\partial \Delta \boldsymbol{\epsilon}}$ ; de la partie déviatorique de la prédiction élastique (4.19a), on comprend que :

$$
\mathbf{s}_{n+1}^* = \mathbf{s}_n + 2\mu \Delta \boldsymbol{\epsilon}^d \quad \text{avec} \quad \boldsymbol{\epsilon}^d = \boldsymbol{\epsilon} - \frac{1}{3} (\operatorname{tr} \boldsymbol{\epsilon}) \mathbf{1} \qquad (4.32)
$$

où $\boldsymbol{\epsilon}^d$ est le déviateur des déformations. En observant que $\boldsymbol{\epsilon}^d = \mathbf{P}^d : \boldsymbol{\epsilon}, \mathbf{P}^d$ étant le projecteur déviatorique, il est aisé de calculer :

$$
\frac{\partial \boldsymbol{s}_{n+1}^*}{\partial \Delta \boldsymbol{\epsilon}} = 2\mu \mathbf{P}^d \qquad (4.33)
$$

2. Calcul de $\frac{\partial \sigma_{eq_{n+1}}^*}{\partial \Delta \boldsymbol{\epsilon}}$ ; considérant la contrainte de von Mises et (4.33), il vient :

$$
\frac{\partial \sigma_{eq_{n+1}}^*}{\partial \Delta \boldsymbol{\epsilon}} = 3\mu \frac{\mathbf{s}_{n+1}^*}{\sigma_{eq_{n+1}}^*} \qquad (4.34)
$$

3. Calcul de $\frac{\partial \Delta p}{\partial \Delta \boldsymbol{\epsilon}}$ ; en dérivant l'équation (4.25) par rapport à $\Delta \boldsymbol{\epsilon}$, il vient :

$$
\frac{\partial \Delta p}{\partial \Delta \boldsymbol{\epsilon}} = \frac{3\mu}{3\mu + \frac{\partial R_{n+1}}{\partial \Delta p}} \frac{\boldsymbol{s}_{n+1}^*}{\sigma_{\text{eq}\,n+1}^*} \quad (4.35)
$$

En tenant compte de (4.33), (4.34) et (4.35), et en les injectant dans (4.31), il vient l'expression
du « correcteur plastique » **E** :

$$
\mathbf{E} = 3\mu(\gamma - \beta)\left(\frac{\mathbf{s}_{n+1}^*}{\sigma_{\text{eq}\,n+1}^*} \otimes \frac{\mathbf{s}_{n+1}^*}{\sigma_{\text{eq}\,n+1}^*}\right) + 2\mu\beta\mathbf{P}^d \quad (4.36)
$$

avec

$$
\gamma = \frac{3\mu}{3\mu + \frac{\partial R_{n+1}}{\partial \Delta p}} \text{ et } \beta = 1 - \frac{R_{n+1}}{\sigma_{\text{eq}\,n+1}^*}
$$
<figure><img src="image_69.png" /></figure>

1: Données : conditions initiales {$ε_n, ε_n^p, (V_k)_n$} et incrément de déformation $Δε$

2: **for** G = 1 **to** N<sub>b</sub> points de Gauss de l'élément **do**

3: Prédiction élastique : $\sigma_{n+1}^* = \sigma_n + \lambda(\text{tr}(\Delta\boldsymbol{\epsilon}))\mathbf{1} + 2\mu\Delta\boldsymbol{\epsilon}$

4: $f_{n+1}^* = \sigma_{eq_{n+1}}^* - hp_n - \sigma_y$

5: **if** $f_{n+1}^* \le 0$ **then**

6: Évolution élastique : actualiser $(\bullet)_{n+1} = (\bullet)_{n+1}^*$ & Sortie

7: **else**

8: Évolution élastoplastique : retour radial

9: Résoudre sur $\Delta p$ l'équation $\sigma_{eq_{n+1}}^* - 3\mu\Delta p - R_{n+1}(p) - \sigma_y = 0$

10: Actualiser les variables :

$$ \Delta\boldsymbol{\epsilon}^p = \frac{3}{2} \frac{\Delta p}{\sigma_{eq_{n+1}}^*} \boldsymbol{s}_{n+1}^* $$

$$ \boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{n+1}^* - 2\mu\Delta\boldsymbol{\epsilon}^p $$

$$ (\bullet)_{n+1} = (\bullet)_n + \Delta(\bullet) $$

11: Calcul du module tangent :

$$ \mathbf{H} = \mathbf{D} - 3\mu(\gamma - \beta)\left(\frac{\mathbf{s}_{n+1}^*}{\sigma_{eq_{n+1}}^*} \otimes \frac{\mathbf{s}_{n+1}^*}{\sigma_{eq_{n+1}}^*}\right) - 2\mu\mathbf{P}^d ; \quad \gamma = \frac{3\mu}{3\mu + \frac{\partial R_{n+1}}{\partial \Delta p}} ; \quad \beta = \frac{3\mu\Delta p}{\sigma_{eq_{n+1}}^*} $$

12: **end if**

13: **end for**

**Algorithmme 4.2 – Algorithme de retour radial**

**Synthèse** L'opérateur tangent est défini en fonction que l'évolution est élastique ou élastoplastique, de sorte que :

$$ \mathbf{H} = \begin{cases} \mathbf{D} & \text{si } f_{n+1}^* \le 0 \text{ évolution élastique} \\ \mathbf{D} - \mathbf{E} & \text{si } f_{n+1}^* > 0 \text{ évolution élastoplastique} \end{cases} \qquad (4.37) $$

L'algorithme 4.2 résume l'implémentation de la méthode du retour radial. Cet algorithme peut faire apparaître une erreur d'intégration temporelle dès que des chargements non proportionnels sont mis en jeu, c'est ce qu'on appelle l'écart à la radialité. En effet, avec la discrétisation implicite utilisée pour cet algorithme, les normales élastique prédite et de la surface de charge finale sont colinéaires. Ceci est vrai quand la normale **n** est constante sur le pas de temps, auquel cas l'intégration est exacte, il y a équivalence entre les lois de comportement continue et discrète. Par conséquent, l'algorithme de retour radial permet une intégration exacte pour un incrément de contrainte radial. Si la normale **n** varie sur le pas de temps $[t_n; t_{n+1}]$ (figure 4.4), on peut définir une erreur induite par cet algorithme et appelée l'écart à la radialité :

$$ \epsilon = \left\| \Delta\boldsymbol{\epsilon} - \sqrt{\frac{3}{2}} \Delta p \mathbf{n}_{n+1} \right\| = \left\| \sqrt{\frac{3}{2}} \int_{t_n}^{t_{n+1}} \dot{p}(\tau) (\mathbf{n}(\tau) - \mathbf{n}_{n+1}) \, d\tau \right\| \quad (4.38) $$

### 4.3.2.2 Méthode des différences centrales généralisées

On généralise l'approche précédente en présentant la méthode des différences centrales généralisées, ou la méthode du point milieu généralisé. Considérons le problème de Cauchy
suivant :

$$
\left.
\begin{aligned}
\dot{X}(t) &= f(X(t)) \\
X(0) &= X_n
\end{aligned}
\right\} \text{dans } [0, T]
\qquad (4.39)
$$

Le principe de base de la méthode des différences centrales généralisées est de supposer
que la vitesse d'une quantité X est constante sur chaque intervalle de temps [t_n; t_{n+1}], et égale
à la valeur à un instant intermédiaire t_{n+θ} ≡ t_θ, c'est-à-dire à t + θΔt (avec 0 ≤ θ ≤ 1) :

$$
\frac{X(t + \Delta t) - X(t)}{\Delta t} = \dot{X}(t + \theta \Delta t) \qquad (4.40)
$$

On écrit alors :

$$
\begin{align}
\Delta X &= X_{n+1} - X_n \nonumber \\
X|_{t_\theta} &= (1-\theta)X_n + \theta X_{n+1} \quad \forall \theta \in [0, 1] \tag{4.41}
\end{align}
$$

On rappelle que si $\theta < 1/2$ la méthode est conditionnellement stable, et si $\theta \ge 1/2$ elle est inconditionnellement stable. Cette méthode permet de retomber sur le retour radial pour $\theta = 1$. Le problème (4.17), discrétisé avec la méthode des différences centrales généralisées s'écrit de la manière suivante :

$$
\begin{align*}
& \sigma_{\text{eq}n+1} - R_{n+1} - \sigma_y = 0 \\
& \Delta \boldsymbol{\epsilon}^p - \Delta p \left. \frac{\mathbf{s}}{\sigma_{\text{eq}}} \right|_{t_\theta} = 0 \tag{4.42} \\
& \Delta \boldsymbol{\sigma} - \mathbf{K} : (\Delta \boldsymbol{\epsilon} - \Delta \boldsymbol{\epsilon}^p) = 0
\end{align*}
$$

où on assure la vérification du critère à la fin de l'intervalle. Ce système peut s'écrire $\mathbb{F}(\mathcal{X}) = 0$, où $\mathbb{F}$ est la fonction résidu à annuler, et $\mathcal{X}$ le vecteur des inconnues tel que :

$$
\mathcal{X} = (\Delta \boldsymbol{\sigma}, \Delta \boldsymbol{\epsilon}^p, \Delta p) \tag{4.43}
$$

On rappelle que l'incrément de déformation $\Delta \boldsymbol{\epsilon}$ est imposé (on pilote l'intégration en déformation). Le système $\mathbb{F}(\mathcal{X}) = 0$ est un système d'équations non linéaires qu'on résoud de façon itérative, par exemple avec une méthode de Newton. Ainsi, en faisant un développement limité de $\mathbb{F}$ au voisinage de $\mathcal{X}^i$ et en cherchant $\mathcal{X}^{i+1}$ tel que $\mathbb{F}(\mathcal{X}^{i+1}) = 0$, l'approximation suivante de la solution est définie par le système linéaire :

$$
\mathbb{F}(\mathcal{X}^{i+1}) \approx \mathbb{F}(\mathcal{X}^i) + \left. \frac{\partial \mathbb{F}(\mathcal{X})}{\partial \mathcal{X}} \right|^{i} \cdot (\mathcal{X}^{i+1} - \mathcal{X}^i) = 0 \quad (4.44)
$$

En faisant le calcul de la matrice jacobienne, on peut observer qu'elle n'est pas symétrique, elle est néanmoins de petite taille. L'algorithme converge quand le résidu $\mathbb{F}(\mathcal{X}^{i+1})$ est inférieur à une valeur donnée par l'utilisateur. On peut remarquer que ce résidu comporte des grandeurs homogènes à des contraintes et à des déformations n'ayant pas le même ordre de grandeur. Il faut donc normaliser ces grandeurs par des grandeurs constantes respectivement homogènes avant le calcul du résidu.

Cet algorithme met aussi en œuvre un état élastique prédit, qui permet d'éviter des calcul
inutiles dans les phases élastiques. La méthode des différences centrales généralisées est
une méthode générale pour l'intégration de relations de comportement, qui a la particularité
d'imposer un chemin pour l'écoulement plastique sur tout l'incrément, ce qui limite forcément
leur taille. Pour obtenir une bonne solution approchée, la direction de l'écoulement plastique
doit varier faiblement sur l'incrément.
<figure><img src="image_70.png" /></figure>

### 4.3.2.3 Le Return Mapping « généralisé »

Un autre algorithme de Return mapping, aussi appelé « cutting-plane algorithm », utilise directement les propriétés de la fonction seuil pour définir le schéma de résolution. Cette méthode passe aussi par un état élastique prédit. Supposons une charge plastique, l'état élastique prédit viole alors le critère de palsticité $f_{n+1}^* > 0$, l'idée principale de cet algorithme est d'effectuer une linéarisation du critère au voisinage de la solution approchée à l'itération i du calcul de l'intégration de la loi de comportement, et de rechercher les corrections permettant d'assurer la nullité du critère à l'instant $t_{n+1}$ :

$$ 
\begin{gathered}
f_{n+1} = \sigma_{eq,n+1} - R_{n+1} - \sigma_y \le 0 \\ 
f_{n+1}^{(i+1)} \approx f_{n+1}^{(i)} + \left. \frac{\partial f}{\partial \boldsymbol{\sigma}} \right|_{n+1}^{(i)} : (\boldsymbol{\sigma}_{n+1}^{(i+1)} - \boldsymbol{\sigma}_{n+1}^{(i)}) + \left. \frac{\partial f}{\partial R} \right|_{n+1}^{(i)} (R_{n+1}^{(i+1)} - R_{n+1}^{(i)}) = 0
\end{gathered} 
\quad (4.45) 
$$

Les algorithmes de type return mapping sont basés sur une prédiction élastique puis une correction plastique, qui exploite la séparation de l'opérateur élastique-plastique. Ceci est en fait une conséquence naturelle du fait que les équations constitutives du modèle peuvent être séparées en parties élastique et plastique[8]. Il s'ensuit alors que pour la partie plastique de ces équations, la configuration spatiale est considérée comme fixée. Cette étape est alors considérée comme une étape de relaxation des contraintes par écoulement plastique. Ceci se traduit par l'hypothèse d'un taux de déformation nul $\dot{\boldsymbol{\epsilon}} = \mathbf{0}$. La direction de recherche retenue pour la correction des variables $\boldsymbol{\sigma}$ et $R$ est la direction associée aux lois d'évolution et aux lois d'état. Les équations de la correction plastique s'écrivent de façon générale :

$$ 
\begin{gathered}
\dot{\boldsymbol{\epsilon}} = \dot{\boldsymbol{\epsilon}}^e + \dot{\boldsymbol{\epsilon}}^P = \mathbf{0} \\ 
\dot{\boldsymbol{\epsilon}}^P = \dot{p} \frac{\partial f}{\partial \boldsymbol{\sigma}} \\ 
\dot{\boldsymbol{\sigma}} = \mathbf{K} : \dot{\boldsymbol{\epsilon}}^e \\ 
\dot{R} = h \dot{p}
\end{gathered} 
\quad (4.46) 
$$

où $\mathbf{K}$ est le tenseur d'élasticité. En linéarisant ces lois et en les combinant, il vient :

$$ 
\begin{gathered}
\sigma_{n+1}^{(i+1)} - \sigma_{n+1}^{(i)} = -\delta p \mathbf{K} : \left. \frac{\partial f}{\partial \boldsymbol{\sigma}} \right|_{n+1}^{(i)} \\ 
R_{n+1}^{(i+1)} - R_{n+1}^{(i)} = h \delta p
\end{gathered} 
\quad (4.47) 
$$

où $\delta p$ est défini comme $\delta p = p_{n+1}^{i+1} - p_{n+1}^{(i)}$. En reportant (4.47) dans (4.45), on déduit l'expression de la correction de la déformation plastique cumulée :

$$ 
\delta p = \frac{f_{n+1}^{(i)}}{\left. \frac{\partial f}{\partial \boldsymbol{\sigma}} \right|_{n+1}^{(i)} : \mathbf{K} : \left. \frac{\partial f}{\partial \boldsymbol{\sigma}} \right|_{n+1}^{(i)} - h \left. \frac{\partial f}{\partial R} \right|_{n+1}^{(i)}}, \quad p_{n+1}^{(i+1)} = p_{n+1}^{(i)} + \delta p \quad (4.48) 
$$

Les itérations se poursuivent jusqu'à vérification du critère à $t_{n+1}$. En considérant un comportement élastique isotrope, l'expression (4.48) se simplifie et donne :

$$ 
\delta p = \frac{f_{n+1}^{(i)}}{h + 3\mu}, \quad p_{n+1}^{(i+1)} = p_{n+1}^{(i)} + \delta p \quad (4.49) 
$$

L'algorithme 4.3 résume la procédure du Return Mapping.

<footnote>[8] SIMO et HUGHES. *Computational Inelasticity*. Springer, 1998.</footnote>