# Une introduction à la plasticité cristalline interactions avec l'environnement

Georges Cailletaud

► To cite this version:

Georges Cailletaud. Une introduction à la plasticité cristalline interactions avec l'environnement. PlastOx 2007, May 2007, Argelès-sur-Mer, France. pp.81-115, 10.1051/ptox/2009008. hal-01497918

HAL Id: hal-01497918

https://hal.science/hal-01497918v1

Submitted on 29 Mar 2017

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Une introduction à la plasticité cristalline interactions avec l'environnement

G. Cailletaud

Centre des Matériaux, CNRS UMR 7633, MINES ParisTech, BP. 87,
91003 Evry Cedex, France

e-mail: Georges.Cailletaud@ensmp.fr

**Résumé.** Cette partie est une introduction aux modèles phénoménologiques de plasticité cristalline. On discute leurs propriétés, et leurs utilisations dans le cadre de modélisations de microstructures. Les versions plastiques et viscoplastiques sont citées. Les applications portent sur des monocristaux ou des polycristaux. Dans ce dernier cas, on a privilégié les modèles éléments finis, qui permettent de véritables calculs de microstructures. L'exposé se conclut par quelques exemples dans lesquels l'environnement jour un rle important.

## 1. INTRODUCTION

Le cadre théorique qui a permis le développement des modèles de monocristal a été construit dans les années 70 [1–6], et les premières applications ont suivi quelques années plus tard [7, 8]. L’histoire de la plasticité cristalline est donc bien connue, elle peut-être consultée dans quelques livres de référence, comme [9, 10]. Le but de cette section n’est pas d’ajouter une nouvelle description de ces approches maintenant relativement classiques, mais de présenter le domaine en suivant une approche alternative, basée sur la thermodynamique des processus irréversibles. Le modèle de monocristal est ainsi vu comme un modèle « multipotentiel ». On illustre ensuite les possibilités de l’approche micromécanique mise en place, dans le cadre de modèles purement mécaniques et dans le cadre de couplages multiphysiques.

L’intérêt de l’approche thermodynamique est de fournir des schémas naturels permettant de choisir les variables critiques des modèles et la forme de leurs lois d’évolution. En l’absence d’endommagement, on suppose que les mécanismes de déformation sont réduits au glissement cristallographique, sur un plan donné et dans une direction donnée. La classe de modèles qui découle de cette hypothèse peut-être utilisée pour la description de monocristaux comme de polycristaux. Le cas du monocristal est illustré en Fig. 1, qui montre le développement de la plasticité dans un monocristal de superalliage base nickel après 1, 2 et 11 cycles : un premier plan est activé après la première traction, à la suite de quoi d’autres plans peuvent-être observés. Le monocristal est orienté pour le glissement simple, si bien qu’il n’y a qu’une ligne visible en surface. Le glissement multiple est également courant dans les monocristaux, mais plus encore dans les polycristaux, en raison des incompatibilités cinématiques entre grains. La figure 2 montre le cas d’un polycristal CFC (Fig. 2a) et celui d’un cristal de microstructure hexagonale (Fig. 2b). Pour le second cas, le nombre de plans de glissement possibles est relativement réduit, en comparaison avec un matériau CFC, si bien que d’autres mécanismes deviennent actifs, comme le maclage ou le dommage intergranulaire.

Les modèles de plasticité cristalline développés au cours des vingt dernières années proviennent de deux communautés de recherche : il y a les modèles « purement phénoménologiques », et ceux qui sont basés sur l’évolution de densités de dislocations. Afin de conserver au texte une taille raisonnable, il a été décidé de mettre en avant l’aspect « opérationnel » des modèles, si bien que la présentation sera centrée sur les modèles phénoménologiques. Cependant un lien sera fait entre les deux approches.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Activation successive des plans de glissement dans une éprouvette monocristalline [14].</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Activité des systèmes de glissement dans des éprouvettes polycristallines : (a) waspaloy, (b) zircaloy, avec présence de dommage intergranulaire.</figcaption></figure>

Le plan de l'exposé est le suivant :

- en partie 2, afin de poser le cadre et les notations, on commence par rappeler brièvement les éléments essentiels de l'approche thermodynamique ; on montre ensuite une classe de modèles de plasticité cristalline, de type phénoménologique ;

- on introduit ensuite (partie 3) les outils numériques qui seront utilisés dans les analyses par éléments finis ;

- les trois dernières parties sont consacrées à des applications, respectivement calcul de structures, calcul de microstructures et calculs multiphysiques, en parties 4 à 6.
## 2. UNE APPROCHE THERMODYNAMIQUE DE LA PLASTICITÉ CRISTALLINE

Le cadre thermodynamique général [11–13] est d'abord rappelé pour le cas des modèles inélastiques. La présentation est réduite au cas des petites déformations, car c'est le domaine d'utilisation typique dans le présent document.

### 2.1 Cadre général

Le premier principe exprime la conservation de l'énergie. Au sein d'un domaine $D$, l'énergie interne $E$, qui est obtenue en intégrant l'énergie interne spécifique $e$ est la somme de la puissance des efforts intérieurs et des échanges de chaleur:

$$ \frac{dE}{dt} = \int_{D} \rho \frac{de}{dt} dV = \mathcal{P}^{(e)} + \dot{Q} \quad (1) $$

avec les définitions suivantes :

- Puissance des efforts intérieurs :

$$ \mathcal{P}^{(e)} = \int_{D} \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} dV \quad (2) $$

- Quantité de chaleur échangée, calculée à partir du taux de chaleur reçue, $\underline{q}$, et de la chaleur volumique, $r$, $\underline{n}$ désignant la normale sortante à la surface $\partial D$ :

$$ \dot{Q} = \int_{D} r dV - \int_{\partial D} \underline{q} \cdot \underline{n} dS = \int_{D} (r - \operatorname{div} \underline{q}) dV \quad (3) $$

Ceci fournit l'expression suivante pour le premier principe :

$$ \rho \frac{de}{dt} = \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} + r - \operatorname{div} \underline{q} \quad (4) $$

Le deuxième principe fournit une borne supérieure de la chaleur reçue par le volume $D$ à une température $T$, et peut être exprimé comme une fonction de l'entropie $S$ ou de l'entropie spécifique $s$ :

$$ \frac{dS}{dt} \geq \int_{D} \frac{r}{T} dV - \int_{\partial D} \frac{\underline{q} \cdot \underline{n}}{T} dS \quad (5) $$

$$ \text{alors :} \quad \int_{D} \left( \rho \frac{ds}{dt} - \frac{r}{T} + \operatorname{div} \left( \frac{\underline{q}}{T} \right) \right) dV \geq 0 \quad (6) $$

En utilisant l'énergie libre de Helmholtz, $\Psi$, telle que $e = \Psi + Ts$, on obtient l'inégalité dite de Clausius-Duhem :

$$ \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} - \rho \frac{d\Psi}{dt} - \rho s \dot{T} - \frac{1}{T} \underline{q} \cdot \underline{\operatorname{grad}}(T) \geq 0 \quad (7) $$

La méthode de l'état local suppose que l'ensemble de l'histoire du matériau peut être représentée par la température, $T$, et par un ensemble de variables d'état $\alpha_i$, qui sont présentes dans $\Psi$, et qui seront utilisées dans l'argument classique de Coleman-Noll. On trouve successivement :

$$ \frac{d\Psi}{dt} = \frac{\partial \Psi}{\partial T} \dot{T} + \frac{\partial \Psi}{\partial \alpha_i} \dot{\alpha}_i \quad (8) $$

$$ s = -\frac{\partial \Psi}{\partial T} \quad (9) $$

$$ \sigma_{ij} \dot{\epsilon}_{ij} - \rho \frac{\partial \Psi}{\partial \alpha_i} \dot{\alpha}_i - \frac{1}{T} \underline{q} \cdot \underline{\operatorname{grad}}(T) \geq 0 \quad (10) $$
La dissipation intrinsèque $\Phi_1$ et la dissipation thermique $\Phi_2$ sont considérées :

$$ \Phi_1 = \sigma_{ij} \dot{\epsilon}_{ij} - \rho \frac{\partial \Psi}{\partial \alpha_i} \dot{\alpha}_i \qquad \Phi_2 = -\frac{1}{T} \underline{q} \cdot \underline{\mathrm{grad}}(T) \qquad (11) $$

On suppose alors que chacune d'entre elles doit être positive. Cette positivité est directement assurée par la loi de Fourier pour ce qui concerne la partie thermique. On obtient ainsi l'équation de la chaleur en présence de déformation mécanique :

$$ \underline{q} = -k(T, \alpha_i) \underline{\mathrm{grad}}(T) \qquad (12) $$

$$ \mathrm{div} \left( k \underline{\mathrm{grad}}(T) \right) = \rho C_{\varepsilon} \dot{T} - r - \sigma_{ij} \dot{\epsilon}_{ij} + \rho \left( \frac{\partial \Psi}{\partial \alpha_i} - T \frac{\partial^2 \Psi}{\partial T \partial \alpha_i} \right) \dot{\alpha}_i \qquad (13) $$

(avec $C_\varepsilon = T \partial s / \partial T$, chaleur massique spécifique à déformation constante). La forme de la partie mécanique de cette équation peut prendre différentes formes, en fonction du type de matériau considéré et des phénomènes que l'on souhaite modéliser.

### 2.1.1 Thermoélasticité

En thermoelasticité, la seule variable interne est la déformation élastique. Dans la mesure où l'élasticité est un processus non dissipatif, $\Phi_1$ doit rester à zéro. Ceci conduit immédiatement à la définition du tenseur de contrainte :

$$ \Phi_1 = \underline{\sigma} : \dot{\underline{\epsilon}}^e - \rho \frac{\partial \Psi}{\partial \underline{\epsilon}^e} : \dot{\underline{\epsilon}}^e = 0 \qquad (14) $$

Les équations de l'élasticité impliquent donc deux variables d'état, $T$ et $\underline{\epsilon}^e$, et deux variables conjuguées, l'entropie $s$ et $\alpha$ tenseur de contrainte $\underline{\sigma}$. $\Psi$ peut être vue comme le potentiel thermodynamique qui caractérise le processus réversible.

### 2.1.2 Processus dissipatifs

On introduit maintenant une série de variables d'écrouissage $\alpha_I$, si bien que le modèle est formulé selon le schéma suivant :

<table><thead><tr><th>Variables d'état</th><th></th><th>Variables conjuguées</th><th></th></tr></thead><tbody><tr><td>température</td><td><i>T</i></td><td><i>s</i> = -<sup>∂<i>Ψ</i></sup>⁄<sub>∂<i>T</i></sub></td><td>entropie</td></tr><tr><td>déformation élastique</td><td><i>ε</i><sup>e</sup></td><td><i>σ</i> = <i>ρ</i> <sup>∂<i>Ψ</i></sup>⁄<sub>∂<i>ε</i><sup>e</sup></sub></td><td>contrainte</td></tr><tr><td>variables d'état</td><td><i>α</i><sub>I</sub></td><td><i>A</i><sub>I</sub> = <i>ρ</i> <sup>∂<i>Ψ</i></sup>⁄<sub>∂<i>α</i><sub>I</sub></sub></td><td>variables d'écrouissage</td></tr></tbody></table>

La dissipation intrinsèque peut alors s'écrire :

$$ \Phi_1 = \underline{\sigma} : \dot{\underline{\epsilon}}^p - A_I \dot{\alpha}_I = Z \dot{z} \qquad (15) $$

avec : $Z = \{\underline{\sigma}, A_I\}$; $z = \{\underline{\epsilon}^p, -\alpha_I\}$.

Un modèle est *standard* [1] si on peut trouver un potentiel $\Omega \equiv \Omega(Z)$ tel que :

$$ \dot{z} = \frac{\partial \Omega}{\partial Z} \qquad (16) $$

Si $\Omega$ est une fonction convexe de $Z$ qui contient l'origine, la dissipation est alors automatiquement positive, puisque :

$$ \phi_1 = Z \frac{\partial \Omega}{\partial Z} \qquad (17) $$
On peut alors définir (au moyen d'une tranformée de Legendre-Fenchel) un potentiel associé, en termes de $\dot{z}$ :

$$ \Omega^*(\dot{z}) = \max_{Z} (Z\dot{z} - \Omega(Z)) \quad (18) $$

En conséquence, $\Omega^*(\dot{z})$ peut être utilisé aussi bien que $\Omega(Z)$ pour caractériser le processus dissipatif. Un modèle dissipatif peut alors être caractérisé complètement par la donnée de deux potentiels, $\Psi$ et $\Omega$ (ou $\Omega^*$). La relation entre les variables d'état et les variables d'écrouissage sont dérivées de $\Psi$, alors que la nature des variables d'écrouissage et leur évolution provient de $\Omega$. A titre d'exemple, un modèle à écrouissage isotrope et cinématique sera obtenu à partir du choix suivant de l'ensemble ($A_I, \alpha_I$) :

<table><thead><tr><th>Type d'écrouissage</th><th>Variable d'état</th><th>Variable conjuguée</th></tr></thead><tbody><tr><td>Ecrouissage isotrope</td><td>r</td><td>R</td></tr><tr><td>Ecrouissage cinématique</td><td>α̃</td><td>X̃</td></tr></tbody></table>

$$ \dot{\tilde{\epsilon}}^p = \frac{\partial \Omega}{\partial \tilde{\sigma}} \qquad \dot{r} = -\frac{\partial \Omega}{\partial R} \qquad \dot{\alpha} = -\frac{\partial \Omega}{\partial X} \quad (19) $$

En supposant que $J$ est la contrainte équivalente de von Mises, telle que $J(\tilde{x}) = ((3/2)\tilde{x} : \tilde{x})^{1/2}$, on obtient simplement un modèle viscoplastique au moyen d'un potentiel $\Omega$, avec les étapes suivantes :

- Définition du seuil de viscoplasticité :

$$ f(\tilde{\sigma}, X, R) = J(\tilde{\sigma} - X) - R \quad (20) $$

- Définition d'un potentiel, fonction de $f$ :

$$ \Omega = \frac{K}{n+1} \left(\frac{f}{K}\right)^{n+1} \quad (21) $$

- Ecoulement viscoplastique :

$$ \dot{\tilde{\epsilon}}^p = \frac{\partial \Omega}{\partial \tilde{\sigma}} = \frac{\partial \Omega}{\partial f} \frac{\partial f}{\partial \tilde{\sigma}} \quad (22) $$

- Intensité de l'écoulement (ici, modèle de Norton avec une contrainte interne) :

$$ \dot{v} = \frac{\partial \Omega}{\partial f} = \left(\frac{f}{K}\right)^n \quad (23) $$

- Direction de l'écoulement (normale à la surface équipotentielle) :

$$ \tilde{n} = \frac{\partial f}{\partial \tilde{\sigma}} \quad (24) $$

- Ecrouissage :

$$ \dot{r} = -\frac{\partial \Omega}{\partial R} \qquad \dot{\alpha} = -\frac{\partial \Omega}{\partial X} \quad (25) $$

On peut passer d'un modèle viscoplastique à un modèle de plasticité indépendante du temps en remplaçant le potentiel viscoplastique par un pseudo-potentiel plastique, comme le montre la Fig. 3, où Ind(f) passe de 0 (pour $f < 0$) à l'infini (pour $f = 0$). La vitesse d'écoulement n'est alors plus déterminée par la dérivée partielle de $\Omega$ en fonction de $f$, mais par la condition de cohérence. La direction de l'écoulement peut également se retrouver en introduisant le principe de Hill, qui suppose que le champ de contrainte réel permet d'atteindre le maximum de la puissance intrinsèque dissipée $\Phi_1 = \tilde{\sigma}^* : \dot{\tilde{\epsilon}}^p$ parmi tous les champs de contrainte admissibles, lorsqu'il n'y a pas d'écrouissage :

$$ \forall \tilde{\sigma}^* \text{ admissible } (\tilde{\sigma} - \tilde{\sigma}^*) : \dot{\tilde{\epsilon}}^p \ge 0 \quad (26) $$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Vue schématique (a) du potentiel viscoplastique, (b) du pseudo-potentiel plastique.</figcaption></figure>

Dans le cas d'un matériau avec écrouissage, l'équation précédente doit être reformulée en :

$$ \Phi_1 = \underline{\sigma}^* : \dot{\underline{\epsilon}}^p - A_I^* \dot{\alpha}_I = \underline{\sigma} : \dot{\underline{\epsilon}}^p - \dot{\Psi}_p = Z \dot{z} \quad \text{maximum} \quad (27) $$

où Z comprend le tenseur de contrainte et les variables d'écrouissage $A_I$, et où z comprend le tenseur de déformation plastique et les variables d'état $(-\alpha_I)$ :

$$ (Z - Z^*) \dot{z} \ge 0 \quad (28) $$

La maximisation de $\Phi_1$ sous la contrainte $f \le 0$ peut être vue comme une extension du principe de Hill. Si on définit $F(Z) = Z\dot{z} - \lambda f$ et qu'on cherche le zéro de $\partial F / \partial Z$

$$ \dot{z} = \lambda \frac{\partial f}{\partial Z} \quad \text{soit :} \quad \dot{\underline{\epsilon}}^p = \lambda \frac{\partial f}{\partial \underline{\sigma}} = \lambda \underline{n} \qquad \dot{\alpha}_I = -\lambda \frac{\partial f}{\partial A_I} \quad (29) $$

$\lambda$ (inconnu pour le moment) joue le même rôle que la vitesse de déformation équivalente en viscoplasticité.

Ce type d'approche offre l'occasion d'introduire deux types de couplage :

- couplage d'état, dans l'énergie libre (à noter la symétrie des interactions) :

$$ \Psi(\alpha_1, \alpha_2) = \frac{1}{2} c_{11} \alpha_1^2 + \frac{1}{2} c_{22} \alpha_2^2 + \frac{1}{2} c_{12} \alpha_1 \alpha_2 \quad (30) $$

$$ \frac{\partial A_1}{\partial \alpha_2} = \frac{\partial A_2}{\partial \alpha_1} = \frac{\partial^2 \phi}{\partial \alpha_1 \partial \alpha_2} \quad (31) $$

- Couplage dissipatif, lorsque $\Omega$ est la somme de plusieurs fonctions potentiel, $\Omega_K$ :

$$ \dot{z} = \sum_K \frac{\partial \Omega_K}{\partial Z} \quad (32) $$

Le lecteur concerné pourra consulter sur le sujet quelques ouvrages de référence, comme [2–4].
## 2.2 Développement d'une classe de modèles de monocristal

### 2.2.1 Cadre général

Le monocristal est vu comme une collection de N systèmes de glissement, définis par leur plan de glissement, de normale $\underline{n}^s$, et leur direction de glissement $\underline{l}^s$. Le tenseur d'orientation s'écrit :

$$ \underline{m}^s = \frac{1}{2}(\underline{n}^s \otimes \underline{l}^s + \underline{l}^s \otimes \underline{n}^s) \qquad (33) $$

En petites perturbations, la cission résolue s'exprime :

$$ \tau^s = \underline{\sigma} : \underline{m}^s \qquad (34) $$

et le tenseur de déformation se décompose en une déformation élastique et une déformation inélastique :

$$ \dot{\underline{\varepsilon}} = \dot{\underline{\varepsilon}}^e + \dot{\underline{\varepsilon}}^p \qquad (35) $$

Le comportement élastique doit prendre en compte l'anisotropie (élasticité cubique ou orthotrope en général), mais reste classiquement une relation biunivoque entre contrainte et déformation élastique. La partie dissipative du modèle doit être exprimée pour chaque système de glissement. Les variables d'écrouissage isotrope et cinématique doivent donc être définies sur chaque système, selon le tableau suivant :

<table><thead><tr><th>Phénomène</th><th>Variable d'état</th><th>Variable associée</th></tr></thead><tbody><tr><td>Elasticité</td><td>ε̃<sup>e</sup></td><td>σ̃</td></tr><tr><td>Ecrouissage isotrope</td><td>ρ<sup>s</sup>, s = 1..N</td><td>r<sup>s</sup>, s = 1..N</td></tr><tr><td>Ecrouissage cinématique</td><td>α<sup>s</sup>, s = 1..N</td><td>x<sup>s</sup>, s = 1..N</td></tr></tbody></table>

Les relations entre variables d'état et variables associées proviennent de l'énergie libre, où l'on découpe la partie élastique et la partie inélastique :

$$ \rho\psi(\underline{\varepsilon}^e, \rho^s, \alpha^s) = \rho\psi^e(\underline{\varepsilon}^e) + \rho\psi^p(\rho^s, \alpha^s) \qquad (36) $$

L'élasticité (processus non dissipatif) est totalement définie par :

$$ \rho\psi^e = \frac{1}{2} \underline{\varepsilon}^e : \underline{\Delta} : \underline{\varepsilon}^e \qquad (37) $$

Les variables d'état $\alpha^s$ et $\rho^s$ sont présentes dans l'énergie libre ; les variables d'écrouissage $x^s$ et $r^s$ sont leurs dérivées partielles :

$$ \rho\psi^p = \frac{1}{2} c \sum_s (\alpha^s)^2 + \frac{1}{2} Q \sum_r \sum_s h_{rs} \rho^r \rho^s \qquad (38) $$

$$ x^r = c\alpha^r; \quad r^r = bQ \sum_s h_{rs} \rho^s \qquad (39) $$

On remarque dans ces équations la *matrice d'interaction*, dont les composantes $h_{rs}$ caractérisent l'auto-écrouissage ($r = s$) et l'écrouissage latent (si $r \neq s$). L'effet « multi-mécanisme » a été étudié en détail dans le passé [5, 6] d'un point de vue théorique.

Pour construire un modèle viscoplastique, il suffit de donner un potentiel viscoplastique ; pour construire un modèle de plasticité indépendante du temps, il faut donner un pseudo-potentiel. Ceux-ci sont construits à partir de l'expression de la fonction-seuil sur chaque système de glissement :

$$ f^s = |\tau^s - x^s| - r^s - \tau_0 \qquad (40) $$
### 2.2.2 Formulation viscoplastique

A titre d'exemple, on choisit ici une fonction puissance pour construire le potentiel viscoplastique, sachant que toute autre fonction croissante serait une bonne candidate :

$$ \Omega = \sum_{r} \Omega_{r}(f^{r}) = \frac{K}{n+1} \sum_{r} \left\langle \frac{f^{r}}{K} \right\rangle^{n+1} \quad (41) $$

L'écoulement viscoplastique s'écrit alors :

$$ \dot{\underline{\epsilon}}^p = \frac{\partial \Omega}{\partial \underline{\sigma}} = \sum_{r} \frac{\partial \Omega_r}{\partial \underline{\sigma}} = \sum_{r} \frac{\partial \Omega_r}{\partial f^r} \frac{\partial f^r}{\partial \underline{\sigma}} = \sum_{r} \dot{v}^r \underline{m}^r \eta^r = \sum_{r} \dot{\gamma}^r \underline{m}^r \quad (42) $$

Dans l'équation précédente, $\eta^r$ désigne $\text{signe}(\tau^r - x^r)$ ; la partie scalaire de l'expression caractérise la vitesse de glissement viscoplastique :

$$ \frac{\partial \Omega}{\partial f^r} = \left\langle \frac{f^r}{K} \right\rangle^n = \dot{v}^r \qquad \dot{\gamma}^r = \dot{v}^r \text{signe}(\tau^r - x^r) = \dot{v}^r \eta^r \quad (43) $$

et la partie tensorielle la direction de glissement, donnée par la règle de normalité :

$$ \frac{\partial f^r}{\partial \underline{\sigma}} = \frac{\partial (|\underline{m}^r : \underline{\sigma} - x^r| - r^r - \tau_0)}{\partial \underline{\sigma}} = \underline{m}^r \eta^r \quad (44) $$

Pour un modèle standard, les variables associées aux écrouissages isotrope et cinématique sont donc imposées par le cadre théorique :

- l'écrouissage cinématique sur le système $s$ est gouverné par le glissement lui-même, $\gamma^s$ :

$$ \dot{\alpha}^s = -\frac{\partial \Omega}{\partial x^s} = -\sum_r \frac{\partial \Omega_r}{\partial x^s} = -\frac{\partial \Omega_s}{\partial f^s} \frac{\partial f^s}{\partial x^s} = \dot{v}^s \eta^s = \dot{\gamma}^s \quad (45) $$

- l'écrouissage isotrope sur le système $s$ est gouverné par le glissement cumulé sur le système, $v^s$ :

$$ \dot{\rho}^s = -\frac{\partial \Omega}{\partial r^s} = -\sum_r \frac{\partial \Omega_r}{\partial r^s} = -\frac{\partial \Omega_s}{\partial f^s} \frac{\partial f^s}{\partial r^s} = \dot{v}^s \quad (46) $$

On définira un modèle non standard, par exemple par :

$$ \dot{\alpha}^s = (\eta^s - d\alpha^s)\dot{v}^s \quad (47) $$

$$ \dot{\rho}^s = (1 - b\rho^s)\dot{v}^s \quad (48) $$

Comme $\underline{\sigma} : \dot{\underline{\epsilon}}^p = \underline{\sigma} : \sum_s \underline{m}^s \dot{\gamma}^s = \sum_s \underline{\sigma} : \underline{m}^s \dot{\gamma}^s = \sum_s \tau^s \dot{\gamma}^s = \sum_s \tau^s \eta^s \dot{v}^s$ la dissipation intrinsèque s'écrit :

$$ \phi_1 = \underline{\sigma} : \dot{\underline{\epsilon}}^p - \sum_s x^s \dot{\alpha}^s - \sum_s r^s \dot{\rho}^s \quad (49) $$

$$ = \sum_s (\tau^s \eta^s \dot{v}^s - x^s (\eta^s - d\alpha^s) \dot{v}^s - r^s (1 - b\rho^s) \dot{v}^s) \quad (50) $$

$$ = \sum_s (f^s + \tau_0 + \frac{d}{c} (x^s)^2 + b r^s \rho^s) \dot{v}^s \quad (51) $$

Ceci met en évidence les contributions suivantes :

- dissipation visqueuse, $\sum_s f^s \dot{v}^s$

- dissipation de friction, $\tau_0 \sum_s \dot{v}^s$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Illustration du fonctionnement (a) d'un modèle viscoplastique, (b) d'un modèle plastique.</figcaption></figure>

– dissipation due à la non-linéarité de l’écrouissage, $(\frac{d}{c} (x^s)^2 + br^s \rho^s) \dot{v}^s$.

### 2.2.3 Formulation plastique

La figure 4 illustre la différence entre formulation viscoplastique et formulation plastique. Comme dans le cas des modèles à monopotentiel classiques, le point de fonctionnement courant dans l'espace des contraintes doit rester sur la surface de charge au lieu de se placer sur une équipotentielle. La vitesse d'écoulement plastique pour chaque système de glissement n'est donc plus définie par la distance au seuil (qui reste nulle), mais doit être déduite de la condition de cohérence. La régularisation qui était assurée par l'effet visqueux n'est plus présente, et la direction de l'écoulement doit être définie par une règle additionnelle, lorsque le point de fonctionnement se trouve à un point anguleux du domaine, c'est-à-dire en situation de glissement multiple. Le point crucial est alors de déterminer le jeu de systèmes qui sont effectivement actifs. Un multiplicateur plastique est défini pour chaque système de glissement. Le jeu de multiplicateurs non nuls doivent être solution du système linéaire formé par la condition de cohérence sur les systèmes actifs, $\dot{f}^s = 0$.

$$ \dot{\underline{\epsilon}}^p = \sum_{s=1}^{N} \dot{\lambda}^s \frac{\partial f^s}{\partial \underline{\sigma}} = \sum_{s=1}^{N} \dot{\lambda}^s \underline{m}^s \qquad (52) $$

Un calcul simple permet d'obtenir le système adéquat, à partir de :

– la décomposition de la déformation :

$$ \dot{\underline{\sigma}} = \underline{\Lambda} : \left( \dot{\underline{\epsilon}} - \sum_{r} \underline{m}^{r} \dot{\gamma}^{r} \right) \qquad (53) $$

– la condition de cohérence appliquée aux systèmes actifs :

$$ f^s = |\underline{m}^s : \underline{\sigma} - x^s| - r^s - \tau_0 = 0 \qquad 0 = \underline{m}^s : \dot{\underline{\sigma}} - \dot{x}^s - \eta^s \dot{r}^s \qquad (54) $$

où l'on utilise la notation $\eta^s = \text{signe}(\tau^s - x^s)$.

L'étape suivante consiste à calculer $\underline{m}^s : \dot{\underline{\sigma}}$ dans l'équation (54), en introduisant $H_{sr}$, comme spécifié plus loin, et en conservant la notation $\dot{v}^s$ pour l'intensité de la vitesse de glissement sur le système r :

$$ \underline{m}^s : \dot{\underline{\sigma}} = \dot{x}^s + \eta^s \dot{r}^s = \sum_{r} H_{sr} \dot{v}^r \qquad (55) $$
On remplace alors $\dot{\boldsymbol{\sigma}}$ par son expression dans l'équation (53):

$$ \boldsymbol{m}^s : \boldsymbol{\Lambda} : \dot{\boldsymbol{\epsilon}} - \sum_r \boldsymbol{m}^s : \boldsymbol{\Lambda} : \boldsymbol{m}^r \eta^r \dot{v}^r = \sum_r H_{sr} \dot{v}^r \quad (56) $$

A vitesse de déformation imposée, N equations sont ainsi définies (N = nombre de systèmes actifs) pour calculer les multiplicateurs plastiques :

$$ \sum_r (\boldsymbol{m}^s : \boldsymbol{\Lambda} : \boldsymbol{m}^r + H_{sr}) \dot{v}^r = \boldsymbol{m}^s : \boldsymbol{\Lambda} : \dot{\boldsymbol{\epsilon}} \quad (57) $$

Plusieurs jeux de systèmes de glissement peuvent produire le même tenseur de vitesse de déformation plastique. Il faut donc ajouter une condition pour sélectionner le bon jeu de systèmes, comme on le discute en section (2.5). Dans le présent modèle, la matrice formée par les coefficients $H_{sr}$ est non symétrique :

$$ \dot{x}^s = c \dot{\gamma}^s - dx^s \dot{v}^s = (c \eta^s - dx^s) \dot{v}^s = \sum_r (c \eta^r - dx^r) \delta_{rs} \dot{v}^r \quad (58) $$

$$ \dot{r}^s = Q \sum_r b h_{sr} \exp(-b v^r) \dot{v}^r \quad (59) $$

$$ H_{sr} = (c\eta^r - dx^r)\delta_{sr} + Q b h_{sr} \exp(-b v^r) \quad (60) $$

Néanmoins, elle devient symétrique s'il n'y a pas d'écrouissage cinématique et un écrouissage isotrope linéaire, $r^s = \sum_r h_{sr} v^r$

$$ H_{sr} = h_{sr} \quad (61) $$

Le produit contracté suivant est utile dans les calculs :

$$ \boldsymbol{m}^s : \boldsymbol{I} = m_{ij}^s \delta_{ij} = \operatorname{trace}(\boldsymbol{m}^s) = m_{ii}^s = n_i^s l_i^s = 0 \quad (62) $$

$$ \boldsymbol{m}^s : \boldsymbol{m}^s = m_{ij}^s m_{ij}^s = \frac{1}{4} (n_i^s l_j^s + l_i^s n_j^s) (n_i^s l_j^s + l_i^s n_j^s) = \frac{1}{2} \quad (63) $$

$$ \boldsymbol{m}^s : \boldsymbol{m}^r = m_{ij}^s m_{ij}^r = \frac{1}{4} (n_i^r l_j^r + l_i^r n_j^r) (n_i^s l_j^s + l_i^s n_j^s). \quad (64) $$

Pour le cas de l'élasticité isotrope, le système (57) peut être simplifié, en effet :

$$ \Lambda_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) \quad (65) $$

$$ \Lambda_{ijkl} m_{kl}^r = \mu m_{ij}^r \quad (66) $$

$$ \boldsymbol{m}^s : \boldsymbol{\Lambda} : \boldsymbol{m}^r = \mu m_{ij}^s m_{ij}^r \quad (67) $$

## 2.3 Surfaces de charge

Les surfaces de charge fournissent une information synthétique sur le comportement des matériaux. Le résultat classiquement obtenu pour un monocristal est illustré ici pour le cas d'un matériau CFC, pour lequel le glissement se produit sur les quatre plans octaédriques {111}, dans la direction (110). La surface de charge est définie par une collection d'hyperplans, et possède une forme élémentaire dans des plans tels que celui de traction-cisaillement (Fig. 5a) ou traction biaxiale (Fig. 5b). Pour chacun des deux cas, le repère du monocristal est confondu avec le repère de
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Surface de charge initiale et en cours d'écrouissage : (a) pour un état de contrainte en traction-cisaillement, (b) pour un état de contrainte biaxial.</figcaption></figure>

chargement. La surface de charge initiale n'est pas confondue avec le critère de Tresca ; celui-ci serait obtenu à l'aide d'un super-monocristal offrant des directions de glissement dans toutes les directions de l'espace. Au fur et à mesure d'un chargement à déformation imposée, le point de fonctionnement dans l'espace des contraintes se déplace vers les coins de la surface de charge. Les calculs effectués ici utilisent une hypothèse de petites perturbations, ce n'est donc pas la rotation des plans de glissement qui produit le phénomène. Celui-ci est le résultat de l'équilibre entre déformation élastique et déformation plastique. En Fig. 5a, on considère deux trajets de chargement, le chemin supérieur correspond à un rapport $\epsilon_{12}^{max}/\epsilon_{11}^{max}$ égal à 0,525 alors que pour le chemin inférieur, ce rapport vaut 0,475. Comme le rapport $\epsilon_{12}^{pmax}/\epsilon_{11}^{pmax}$ reste toujours égal à 0,5, il y a une déviation régulière au niveau de la vitesse de contrainte, ce qui explique la forme du trajet dans le plan des contraintes. Le même argument est valide pour le chemin de la Fig. 5b. On note en Fig. 5a que de nouveaux systèmes commencent à s'activer au cours de la déformation. Ce phénomène est réglé par l'écrouissage latent.

## 2.4 Identification pour des chargements de traction et de traction-torsion

Pour une orientation cristalline donnée, le présent modèle peut être comparé à un modèle macroscopique. Son identification peut donc être effectuée par ajustement sur des courbes de comportement, ou encore directement à partir des coefficients d'un modèle macroscopique. Des exemples d'une telle approche peuvent être trouvés dans [7, 8] pour le cas d'un superalliage monocristallin. Si on considère une position de glissement multiple, avec N systèmes de glissement equivalents, et un facteur de Schmid m, on peut transformer l'expression générale donnant $\dot{\epsilon}^p$ :

$$ \dot{\epsilon}^p = \sum_s m_s \dot{v}^s = \sum_s m_s \left\langle \frac{|\tau^s - x^s| - r^s}{k} \right\rangle^n \quad (68) $$

$$ = Nm \left\langle \frac{m(\sigma - x) - r}{k} \right\rangle^n \quad (69) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Simulation of the tensile behaviour for various crystallographic orientations.</figcaption></figure>

qui peut alors être comparée à un modèle macroscopique :

$$ \dot{\epsilon}^p = \left\langle \frac{(\sigma - X) - R}{K} \right\rangle^n \quad \text{with} \quad K = \frac{k}{m} \frac{1}{mN}^{1/n} \quad X = \frac{x}{m} \quad R = \frac{r}{m} \qquad (70) $$

L'exemple de la Fig. 6 montre ainsi les résultats obtenus lors de la simulation d'essais de traction dans 4 directions cristallographiques différentes. On a choisi une valeur négative du paramètre matériau $Q$, afin de montrer les possibilités du modèle.

## 2.5 Sélection des systèmes de glissement

Comme cela a été mentionné précédemment (paragraphe 2.2.3), il n'y a pas unicité du jeu de systèmes de glissement qui conduit à un tenseur vitesse de déformation donné. La condition à remplir pour qu'un ensemble de systèmes soit admissible peut être précisée par approche variationnelle, en utilisant une approche en contrainte ou en vitesse de déformation.

Une façon de poser le problème est de rechercher un jeu de vitesses de glissement admissibles $\dot{v}^s$ qui assure une vitesse de déformation plastique imposée $\dot{\epsilon}^p$ (en supposant $\epsilon \equiv \epsilon^p$). Ceci peut être fait en minimisant la puissance des efforts internes de l'élément de volume [9]:

$$ \mathcal{P}_i = \sum_s \tau_c^s \dot{v}^s \qquad (71) $$

les vitesses de glissement étant soumises à la contrainte :

$$ \mathbf{g}(\dot{v}^s) = \dot{\epsilon}^p - \sum_s \mathbf{m}^s \dot{v}^s \eta^s = 0 \qquad (72) $$

Pour cela, on définit le lagrangien $\mathcal{F}_i$, et on cherche le point-selle :

$$ \mathcal{F}_i(\dot{v}^s, \boldsymbol{\lambda}) = \mathcal{P}_i + \boldsymbol{\lambda} : \mathbf{g} \qquad (73) $$

$$ \frac{\partial \mathcal{F}_i}{\partial \boldsymbol{\lambda}} = \dot{\epsilon}^p - \sum_s \mathbf{m}^s \dot{v}^s = 0 \qquad \frac{\partial \mathcal{F}_i}{\partial \dot{v}^s} = \tau_c^s - \boldsymbol{\lambda} : \mathbf{m}^s \eta^s = 0 \qquad (74) $$
Le tenseur $\lambda$ n'est rien d'autre que le tenseur de contrainte. Pour trouver le jeu de vitesses de glissement qui minimise la puissance interne, il est donc nécessaire de trouver le tenseur de contraintes qui respecte le critère de plasticité, c'est-à-dire qui permet de construire une cission résolue qui atteint $\tau_c^s$ sur les systèmes actifs et qui lui est inférieure sur les systèmes inactifs.

Le problème peut également être posé en choisissant la contrainte comme inconnue, et formulée comme la recherche de l'état de contrainte qui maximise la puissance des efforts extérieurs [10] :

$$ \mathcal{P}_e = \underline{\sigma} : \dot{\underline{\epsilon}}^p \qquad (75) $$

sous la contrainte :

$$ g^s(\underline{\sigma}) = \tau_c^s - \underline{\sigma} : \underline{m}^s \ge 0 \qquad (76) $$

Le lagrangien est alors :

$$ \mathcal{F}_e(\underline{\sigma}, \dot{\mu}^s) = \mathcal{P}_e + \sum_s \dot{\mu}^s g^s \qquad (77) $$

$$ \frac{\partial \mathcal{F}_e}{\partial \underline{\sigma}} = \dot{\underline{\epsilon}}^p - \sum_s \underline{m}^s \dot{\mu}^s = 0 \qquad \frac{\partial \mathcal{F}_e}{\partial \dot{\mu}^s} = \tau_c^s - \underline{\sigma} : \underline{m}^s = 0 \qquad (78) $$

Les multiplicateurs $\dot{\mu}^s$ sont les glissements plastiques. On obtiendra le tenseur de contrainte voulu en trouvant un ensemble de vitesses de glissement qui sont nulles sur les systèmes inactifs, et positives sur les systèmes actifs.

Les deux approches se soldent par la double inégalité [11] :

$$ \underline{\sigma}^* : \dot{\underline{\epsilon}}^p \le \underline{\sigma} : \dot{\underline{\epsilon}}^p = \sum_s \tau_c^s \dot{v}^s \le \sum_s \tau_c^s \dot{v}^{s'} \qquad (79) $$

En résumé :

- parmi tous les états de contrainte admissibles conjugués à la vitesse de déformation plastique réelle, l'état de contrainte réel maximise la puissance plastique. Pour les points réguliers de la surface de charge, ceci indique que la vitesse de déformation plastique est portée par la normale à la surface de charge. Cependant, ceci ne procure pas une information suffisante pour les coins de la surface, en condition de glissement multiple ;

- parmi tous les jeux de vitesses de glissement, le jeu réel minimise la puissance des efforts internes. La dissipation calculée à l'aide l'approche thermodynamique pour un modèle simple comme celui-ci n'a que deux termes :

$$ \Phi_1 = \underline{\sigma} : \dot{\underline{\epsilon}}^p - \sum_s r^s \dot{v}^s \qquad (80) $$

où les variables $r^s$, définies dans l'équation (59), représentent l'augmentation de la cission résolue critique. Comme le montre l'équation (40) (avec $x^s = 0$, puisqu'il n'y a pas d'écrouissage cinématique ici), si tous les systèmes sont équivalents, l'équation se réduit à :

$$ \Phi_1 = \sum_s \tau_0^s \dot{v}^s = \tau_0 \sum_s \dot{v}^s. \qquad (81) $$

## 2.6 Quelques autres modèles de plasticité cristalline

Le modèle qui a été considéré jusqu'à présent représente la branche des modèles de monocristal inspirés des formulations macroscopiques, comme par exemple [12]. Une autre classe de modèles est représentée par les modèles basés sur les densités de dislocations, comme par exemple [13]. D'un autre côté, la formulation de la plasticité cristalline est une solution qui a été utilisée de façon pragmatique pour
représenter d'autres phénomènes qui ont lieu au niveau du réseau cristallin, comme le maclage. Ce paragraphe présente rapidement chacune de ces possibilités.

### 2.6.1 Modèles phénoménologiques et modèles basés sur les densités de dislocations

On parlera de « modèles phénoménologiques » lorsque les variables d'écrouissage sont de type glissement, et de « modèles basés sur les densités de dislocations » lorsque les mécanismes de plasticité liés aux densités de dislocations sont explicitement introduits. Le premier modèle phénoménologique est dû à Taylor [9]. Celui-ci suppose que la cission résolue critique sur les systèmes de glissement dépend du glissement sur l'ensemble des systèmes, par l'intermédiaire d'une matrice d'interaction dont tous les termes sont éaux à 1 : auto-écrouissage et écrouissage latent sont exactement équivalents. Des mesures expérimentales de l'écrouissage latent ont permis d'obtenir de meilleures estimations de la matrice d'écrouissage [14–17]. Celle-ci est maintenant estimée au moyen de calculs de dynamique des dislocations [18]. Le modèle phénoménologique le plus répandu s'exprime dans un cadre viscoplastique, et ne possède pas de seuil [19], si bien que tous les systèmes sont actifs. Il utilise un écrouissage isotrope multiplicatif, qui dépend du glissement cumulé :

$$ \dot{\gamma}^s = \dot{\gamma}_0 \left( \frac{\tau^s}{\tau_c^s} \right)^{1/m} \qquad (82) $$

$$ \tau_c^s = \sum_{r=1}^{N} h_{sr} \dot{\gamma}^r \qquad (83) $$

En fait, les modèles phénoménologiques présentés ici et les modèles basés sur les dislocations ont la même nature : leur variables représentent de façon globale les densités locales de dislocations. Pour illustrer cette affirmation, il suffit de considérer les équations (84)–(86), qui correspondent au modèle [20], dans lesquelles $b$ est le vecteur de Burgers, et $\rho^s$ la densité de dislocations sur le système $s$ :

$$ \dot{\gamma}^s = \dot{\gamma}_0 \left( \frac{\tau^s}{\tau_\mu^s} \right)^{1/m} \qquad (84) $$

$$ \tau_\mu^s = \alpha \mu b \left( \sum_{r=1}^{N} h_{sr} \rho^r \right)^{1/2} \qquad (85) $$

$$ \dot{\rho}^s = \frac{1}{b} \left( \frac{\left( \sum_{r=1}^{N} a^{sr} \rho^r \right)^{1/2}}{K} - 2y_c \rho^s \right) \dot{\gamma}^s \qquad (86) $$

En choisissant $r^s = b^2 \rho^s$ et $y_c = bk_c$, les deux dernières équations peut être réexprimées :

$$ \tau_\mu^s = \alpha \mu \left( \sum_{r=1}^{N} h_{sr} r^r \right)^{1/2} \qquad (87) $$

$$ \dot{r}^s = \left( \frac{\left( \sum_{r=1}^{N} a^{sr} r^r \right)^{1/2}}{K} - 2k_c r^s \right) \dot{\gamma}^s \qquad (88) $$

Les variables $r^s$ qui sont définies ici jouent exactement le rôle de $\rho^s$ dans les equations (38) et (48).
### 2.6.2 Maclage et transformation de phases

Le maclage et la transformation martensitique peuvent être représentées dans un reprère local $\underline{n}^k$, $\underline{t}^k$, $\underline{u}^k$ par le tenseur gradient de transformation $F^k$

$$F = I + \underline{t}^k \otimes \underline{n}^k \qquad (89)$$

Le tenseur de déformation libre est alors $\frac{1}{2}(F^T F - I)$ en grandes transformations, ou encore $\frac{1}{2}(F^T + F) - I$ dans un cadre de petites perturbations. Dans ce dernier cas, la transformation du réseau cristallin produit par le variant $k$ peut être vue comme une déformation paramétrique :

$$\underline{\epsilon}^{tp} := \begin{pmatrix} \epsilon & \gamma & 0 \\ \gamma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad (90)$$

Différentes théories sont disponibles pour déterminer les vecteurs $\underline{n}$ et $\underline{t}$, qui sont fonctions de la cristallographie du matériau (par exemple la théorie de Wechsler-Lieberman-Read, qui fournit 24 variants).

## 3. LES OUTILS DU CALCUL DE MICROSTRUCTURE

Avant de montrer quelques exemples d'utilisation, on rappelle ici quelques méthodes numériques attachées aux calculs avec des modèles de type plasticité cristalline, ainsi que les outils qui sont utiles pour effectuer des calculs prenant en compte l'action de l'environnement. Après un bref énoncé des algorithmes d'intégration particuliers pour le calcul en implicite, on considère les éléments d'interface, et également les méthodes numériques pour les problèmes couplés.

### 3.1 Algorithme d'intégration numérique

Si l'on compare les modèles de monocristaux aux équations des modèles inélastiques classiques, les modèles de monocristal doivent résoudre le problème du choix du jeu de systèmes actifs. C'est pourquoi, après les solutions viscoplastiques classiques dans lesquelles tous les systèmes sont actifs [21, 22], de nombreux auteurs ont proposé des solutions originales [23–28]. On utilise ici une approche implicite de l'intégration, dont l'intérêt est de fournir « gratuitement » la matrice jacobienne nécessaire à la construction de la matrice tangente [4]. En effet, on applique la méthode de Newton pour résoudre le système non linéaire qui provient de l'application de la $\theta$-méthode. Pour des raisons de simplicité, on reporte ici les équations obtenues en petites déformations et écrouissage isotrope linéaire. Les variables d'intégration sont alors les incréments de déformation élastique, $\Delta\epsilon^e$ et les incréments des glissements (visco)plastiques sur tous les systèmes, $\Delta v^s$, qui sont positifs ou nuls. Dans le cas de la plasticité indépendante du temps, ces variables sont déduites de deux types d'équations :

- la décomposition de la déformation :

$$\Delta\underline{\epsilon} = \Delta\underline{\epsilon}^e + \sum_r \underline{m}^r \Delta v^r \qquad (91)$$

- l'expression du critère pour les systèmes de glissement actifs :

$$f^s = |\tau^s| - r^s = 0 \qquad (92)$$

$$\text{with } \tau^s = \underline{m}^s : \underline{\Lambda} : \underline{\epsilon}^e \qquad r^s = \sum_r h_{sr} v^r \text{labeleq : hardx} \qquad (93)$$
Le système résultant est donc formé par les équations (94) et (95):

$$ \mathcal{F}_e = -\Delta\underline{\boldsymbol{\varepsilon}} + \Delta\underline{\boldsymbol{\varepsilon}}^e + \sum_r \underline{\boldsymbol{m}}^r \Delta v^r \quad (94) $$

$$ \mathcal{F}_p^s = \left| \underline{\boldsymbol{m}}^s : \underline{\boldsymbol{\Lambda}} : (\underline{\boldsymbol{\varepsilon}}^e + \Delta\underline{\boldsymbol{\varepsilon}}^e) \right| - \tau_0 - \sum_r h_{sr}(v^r + \Delta v^r) \quad (95) $$

Si le modèle est viscoplastique au lieu d'être plastique, les équations (95) doivent être remplacées par (96), l'équation qui provient de la décomposition de la déformation restant inchangée :

$$ \mathcal{F}_p^s = \left| \underline{\boldsymbol{m}}^s : \underline{\boldsymbol{\Lambda}} : (\underline{\boldsymbol{\varepsilon}}^e + \Delta\underline{\boldsymbol{\varepsilon}}^e) \right| - \tau_0 - \sum_r h_{sr}(v^r + \Delta v^r) - K \left( \frac{\Delta v^s}{\Delta t} \right)^{1/n} \quad (96) $$

Cet algorithme est implanté dans la librairie Zmat du code ZèBuLoN [29] qui a été utilisée pour les simulations réalisées dans ce chapitre.

## 3.2 Eléments d'interface

Le concept, qui est relativement ancien [30, 31], a d'abord été considéré comme une simple solution pour prendre en compte la zone plastique au voisinage de la pointe d'une fissure, afin de dépasser les solutions classiques issues de la théorie de l'élasticité. L'idée a été appliquée plus récemment pour construire un élément fini [32], dit « de zone cohésive ». L'approche est intermédiaire entre la mécanique de la rupture et la mécanique de l'endommagement. Le formalisme peut être dérivé de l'étude d'une fine couche dont on fait tendre l'épaisseur vers zéro [33–35]. L'élément ne possède pas d'épaisseur physique, les nœuds de part et d'autre de l'interface sont confondus deux à deux. Il se décline en 2D (4 nœuds pour les éléments linéaires, 6 nœuds pour les éléments quadratiques) et en 3D (par exemple 6 nœuds pour une face triangulaire linéaire). Chaque nœud porte ses propres déplacements comme inconnues, et les lois de comportement, sont écrites directement en termes de force-déplacement : on introduit donc naturellement l'ouverture, ou déplacement normal au plan de l'élément, et un (2D) ou deux (3D) déplacement tangents, ainsi que les forces nodales conjuguées. Dans les modèles les plus avancés, la courbe de comportement fait intervenir des variables qui prennent en compte l'endommagement et l'effet Bauschinger [34, 36–38]. Comme le précise une revue classique [39], cette classe de modèle peut s'appliquer à n'importe quel type d'interface, pour laquelle l'évolution de l'endommagement va produire une transition entre l'état initial, dans lequel l'élément est élastique, et l'état rompu, où il ne peut plus supporter aucun effort (D = 1). Une formulation de base fait intervenir le comportement sous la forme suivante :

- Elasticité et écoulement plastique :

$$ \begin{align*} \text{saut de déplacement} \qquad & \underline{\boldsymbol{u}} = \underline{\boldsymbol{u}}_e + \underline{\boldsymbol{u}}^p \\ \text{comportement élastique} \qquad & \underline{\boldsymbol{T}} = \underline{\boldsymbol{C}}_e \cdot \underline{\boldsymbol{u}}_e \\ \text{écoulement} \qquad & \dot{\underline{\boldsymbol{u}}}_p = \dot{u}_{eq} \frac{\partial \Phi}{\partial T} \quad \dot{u}_{eq} = \dot{u}_0 \langle \Phi \sigma_0 \rangle^n \end{align*} $$

- Lois d'écoulement :

$$ \begin{align*} \text{cnotrainte équivalente} \qquad & T_{eq} = \frac{1}{F(D)} (\langle T_n \rangle^2 + \beta(T_{t1}^2 + T_{t2}^2))^{0.5} \\ \text{fonction de charge} \qquad & \Phi = T_{eq} - T_0 \end{align*} $$

- Description de l'endommagement :

$$ d_{eq} = u_{eq} / u_{eq}^R \qquad F(d_{eq}) = (1 - d_{eq})^\alpha $$
L'élément fait intervenir les variables suivantes :

- degrés de liberté : $U_{top}$, $U_{bottom}$, ($C_{top}$, $C_{bottom}$)

- saut de déplacement : $\underline{\delta} = N_u(U_{top} - U_{bottom})$

- forces in the element :

$$F_{top}^{u} = \int_{V} N_{u}^{T} \underline{T} dV \qquad F_{top}^{c} = \int_{V} \left( G_{d}^{T} \cdot \underline{j}_{i} + N_{d}(\Delta c_{i} / \Delta t) \right) dV$$

### 3.3 Couplages multiphysiques

La modélisation numérique des problèmes couplés s'effectue en mettant en œuvre dans la procédure numérique un couplage fort, ou un couplage faible. On nomme couplage fort la méthode de résolution employée lorsque l'on résout tous les problèmes en même temps pour un incrément de temps donné, et couplage faible celle qui résulte d'une résolution séquentielle de l'un et l'autre des problèmes. On peut renforcer cette deuxième méthode au moyen d'un point fixe. La première solution nécessite une intervention lourde dans le code de calcul, car il faut développer un élément spécial, à partir d'une formulation variationnelle qui inclut chaque phénomène, et qui aboutit à une écriture dans laquelle on ajoute les degrés de liberté caractérisant chaque problème physique. La matrice de « rigidité » comporte les blocs correspondant à chaque type de problème, et également des blocs d'interaction, qui caractérisent la façon dont les phénomènes sont couplés. Ainsi, un couplage mécanique-diffusion sera caractérisé par un vecteur d'inconnues nodales comportant les déplacements et la concentration. Les blocs de la matrice seront un bloc purement mécanique $[K]_{uu}$, un bloc diffusion, $[K]_{cc}$, et deux blocs d'interaction, $[K]_{uc}$ et $[K]_{cu}$. Ce mode de résolution est présent dans certains codes de calcul commerciaux. Il peut devenir difficile lorsque le conditionnement des matrices se détériore, et nécessite alors des méthodes de résolution adaptées.

La méthode de couplage faible offre une bien plus grande facilité de mise en œuvre, puisqu'il n'est pas nécessaire d'intervenir sur les solveurs. Il n'y a qu'à réaliser l'infrastructure d'échange de données entre différents processus de calcul au niveau de chaque pas. On peut à volonté coupler deux ou plusieurs codes, éventuellement sur des maillages différents, et avec des tailles de pas de calcul différentes selon le problème considéré. Les difficultés rencontrées sont liées à la reprojection des résultats d'un maillage sur l'autre, lorsque ceux-ci sont effectivement différents, et également l'instabilité de la solution si le point fixe n'est pas effectué de façon suffisamment rigoureuse. En fait, cette méthode peut s'avérer recommandable lorsque la physique des phénomènes étudiés implique de toutes façons de faire des petits incréments de temps. La figure 7 montre comment effectuer dans ce cadre un calcul couplé mécanique-diffusion, le couplage étant réalisé par un algorithme de point fixe entre l'instant $t$ et l'instant $t + \Delta t$.

### 4. CALCUL DE STRUCTURES MONOCRISTALLINES

On a à utiliser les modèles de monocristal tels quels dans le cadre de calculs de structures pour des pièces industrielles comme les aubes de turbine, ou encore pour des éprouvettes de laboratoire, dans la mesure où l'aspect monocristallin conduit à des résultats non classiques par rapport aux références classiques sur matériau isotrope.

La figure 8 montre le résultat final d'une chaîne de calculs 3D qui part du calcul thermique, prenant en compte les apports liés aux gaz chauds et le flux de refroidissement dans les canaux internes, puis le calcul mécanique cyclique prenant en compte les camps thermiques, la force centrifuge et les forces aérodynamiques, et enfin le calcul de durée de vie. Les modèles de durées de vie introduisent de façon macroscopique le couplage mécanique-environnement [40] dans le modèle d'endommagement. Il y a en effet trois contributions dans le modèle de durée de vie, une qui est liée à la fatigue, la deuxième au fluage, la troisième à l'oxydation. On pourra consulter [41] pour plus de détails sur la procédure.
<figure><img src="image_8.png" /></figure>

* M et D sont respectivement les variables des problèmes de mécanique et de diffusion

* `mat_param_M` et `mat_param_D` sont les paramètres matériau correspondants

* Les variables diffusion sont vues comme des variables externes pilotées pour le problème mécanique

* Les variables mécaniques sont vues comme des variables externes pilotées pour le problème diffusion

**Figure 7.** Illustration de la méthode de point fixe mise en place pour la résolution d'un problème couplé mécanique-diffusion.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) Prévision de durée de vie sur une aube de turbine subissant un chargement thermomécanique.</figcaption></figure>

Le deuxième exemple rappelle que le caractère monocristallin introduit des hétérogénéités sensibles dans des éprouvettes soumises à des chargements pourtant élémentaires. Une étude détaillée accompagnée de résultats expérimentaux se trouve dans [42]. La figure 9 montre ainsi des cartes de déformation cumulée sur un tube mince ayant subi trois types de chargement :

- une traction simple, l'axe du tube étant confondu avec la direction cristallographique {123} : la face inférieure étant bloquée, on impose un déplacement axial sur la face supérieure, qui ne peut donc pas
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Exemples d'hétérogénéités produites sur un tube monocristallin en traction ou en torsion : (a) traction simple sur un tube dont l'axe est confindu avec la direction cristallographique {123} ; (b) torsion, axe du tube selon la direction {001} ; torsion, axe du tube selon la direction {111}.</figcaption></figure>

tourner, mais qui se déporte hors de l'axe. La déformation plastique est très hétérogène, on retrouve
là une propriété bien connu du monocristal, qui doit faire considérer avec circonspection les essais
réalisés sur des machines trop ridiges et avec des éprouvettes dont la partie utile est courte ;

- une torsion simple, l'axe du tube étant confondu avec la direction cristallographique {001}. Ici, l'éprouvette reste parfaitement dans l'axe, mais, comme le cisaillement macroscopique n'a pas toujours la même direction par rapport aux axes cristallographiques secondaires, on trouve le long de la circonférence une alternance de zones « dures » (en face des directions {100} et {010}) et « molles » (en face des directions {110} et {110});

- une torsion simple, l'axe du tube étant confondu avec la direction cristallographique {111}. Comme précédemment, l'éprouvette reste dans l'axe, mais on retrouve maintenant six zones de déformation préférentielle au lieu de quatre, en raison de la nouvelle symétrie du problème.

Le dernier exemple concerne le développement de la zone plastique en pointe de fissure : ce type de configuration est largement utilisé dans l'étude des modèles locaux d'interaction plasticité-environnement. La première étude du problème [44] a été effectuée de façon semi-analytique sur un matériau parfaitement plastique. La Fig. 10 montre une vérification de ces résultats, au travers d'expérience sur un superalliage base nickel, et par des calculs éléments finis, sur un matériau écrouissable. On observe bien la présence de zones angulaires dans lesquelles la plasticité est totalement absente, même à proximité immédiate de la pointe. L'exposé détaillé de cette étude se retrouve dans [43].

5. DU MONOCRISTAL AU POLYCRISTAL

5.1 Stratégie de transition d'échelles

Dans la littérature, les contraintes internes sont classées en niveaux I, II et III. Le niveau III correspond
à l'échelle microscopique, il capture les fluctuations à l'intérieur des grains, entre l'intérieur et les
zones proches des joints. Le niveau II représente le résultat en moyenne dans une phase (définie par
sa chimie, ou la cristallographie) : les tenseurs résultants en termes de contraintes ou de déformations
ne sont donc plus définis dans chaque grain de façon locale, ils ne prennent pas en compte les relations
d'équilibre locales, ils ne prennent pas en compte le voisinage. Le niveau I se réfère aux contraintes et
aux déformations macroscopiques. Il y a des modélisations numériques attachées à chaque échelle. Les
modèles macroscopiques s'appliquent bien sûr au niveau I. Le niveau II est représenté par les modèles à
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Zone plastique en pointe de fissure. (a) Résultat experimental ; (b) Simulation numerique [43].</figcaption></figure>

champs moyens. Le calcul de microstructures par éléments finis avec des modèles de plasticité cristalline s'appliquent pour évaluer les champs du niveau III.

Pour des raisons d'espace disponible, on ne développera pas les modèles à champs moyens ; la suite du texte va plutôt privilégier la plasticité cristalline par éléments finis, car c'est cette voie qui permet le meilleur dialogue avec les études microstructurales expérimentales classiques. On commente cependant ici la façon dont un assemblage de grains va pouvoir naturellement reproduire le comportement global d'un matériau polycristallin : on se contente de montrer les surfaces de charge, qui représentent en quelque sorte le fondement des modèles, en l'absence d'écrouissage et de considération sur la nature du comportement (plastique ou viscoplastique).

On considère successivement le cas d'un matériau à solidification dirigée (DS) et un polycristal comportant des orientations aléatoires. Afin de simplifier le calcul, on suppose que l'élasticité est uniforme dans tous les grains, en conséquence de quoi il n'y a pas de contraintes internes pendant la phase élastique, donc chaque grain est soumis à la même contrainte. On construit alors simplement les surfaces de charge en considérant tous les systèmes présents dans la microstructure. Pour le matériau DS, tous les grains ont le même axe (001), et une orientation aléatoire dans le plan perpendiculaire à (001). L'axe $x_3$ du repère dans lequel est défini le chargement est confondu avec l'axe cristallographique (001). Le modèle résultant pour un tel agrégat présente encore des caractéristiques de monocristal dans la plan $\sigma_{33}-\sigma_{31}$ (Fig. 11c), en raison du nombre réduit de systèmes de glissement qui sont disponibles pour le cisaillement. D'un autre côté, la forme du critère dans le plan $\sigma_{11}-\sigma_{12}$ (Fig. 11b) est celle d'un critère macroscopique. La figure a été tracée avec trois orientations uniquement, dont les axes (100) font respectivement des angles de 0°, 30°, 60° avec $x_1$. Les références à Tresca et Mises sont tracées en supposant que la cission résolue critique de la loi de Schmid pour chaque monoristal est prise comme cisaillement critique pour le critère de Tresca (différence entre les deux contraintes principales extrêmes), et comme cisaillement octaédral pour le critère de von Mises (la racine carrée de la moitié de la somme des carrés des trois différences entre les différentes contraintes normales principales). Une caractérisation complète du comportement de ce type de matériau a été présentée par ailleurs [45].

Si au lieu de considérer un petit nombre de grains avec des orientations particulières on considère un grand nombre de grains avec des orientations aléatoires (Fig. 11d-f), le critère résultant devient isotrope, et on retrouve, exactement le critère de Tresca (exemple avec 1000 grains).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Surface de charge initiale pour un polycristal. Matériau DS : (a) géométrie de l'agrégat, (b) surface de charge dans le plan $σ_{11}-σ_{12}$, (c) surface de charge dans le plan $σ_{33}-σ_{31}$. Agrégat polycristallin isotrope : (d) géométrie de l'agrégat, (e) surface de charge dans le plan $σ_{11}-σ_{12}$ plane, (f) surface de charge dans le plan $σ_{33}-σ_{31}$ plane.</figcaption></figure>

## 5.2 Le calcul de microstructures

Il y a trois types d'approches en plasticité cristalline par éléments finis, dont les buts diffèrent. En mettant en place les premiers calculs, les auteurs ont voulu faire intervenir le caractère cristallin dans le but de retrouver un comportement macroscopique homogénéisé découlant directement des propriétés locales. L'exigence n'est pas très grande, dans la mesure où on a montré que le comportement global peut s'atteindre de façon fiable avec une représentation très fruste de la microstructure, et un nombre peu élevé de grains. Les premières tentatives remontent aux années 70 [46], mais il a fallu attendre les nouvelles générations de machines pour effectuer des calculs vraiment significatifs, par exemple pour évaluer les évolutions de textures, en 2D ( [47]) ou en 3D ( [48]) : la stratégie est dans ce cas « un élément fini = un grain ». On s'est néanmoins rapidement aperçu que les résultats obtenus ne se démarquent pas notablement du modèle de Taylor (déformation locale dans chaque grain égale à la déformation globale), car un élément fini n'a pas assez de degrés de liberté pour représenter la complexité des champs locaux. Des maillages de plus en plus complexes ont donc été considérés, en introduisant de plus en plus d'éléments pour représenter chaque phase, dans un motif devenant périodique, comme par exemple dans [49] pour les aciers austéno-ferritiques, ou d'autres auteurs qui se sont intéressés à l'intérieur du grain, ou aux zones proches des joints de grains [50–52].

Un autre axe de recherche consiste à s'intéresser au comportement local, en travaillant sur des éprouvettes qui possèdent un nombre réduit de grains. C'est par exemple ce que l'on trouve dans [53, 54], qui effectuent des comparaisons entre champs locaux simulés et expérimentaux. Les études traditionnelles concluent en général à un bon agrément entre simulation et expérience, surtout pour les grains de grande taille. Les joints de grains restent des zones mal caractérisées, en particulier en raison de l'apparition de glissements secondaires. Le comportement d'un grain dans son environnement est différent de celui du monocristal massif du même matériau. Les grains qui sont en sous-couche affectent considérablement le comportement de la surface. Ce point a été étudié en détail de façon
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Carte de déformation plastique équivalente sur un agrégat en cuivre OFHC pour une déformation macroscopique de 5% (valeurs de déformation locale dans le domaine 0–0.214) [57].</figcaption></figure>

numérique dans [55, 56], et aussi dans [57], d'où est extraite la planche de la Fig. 12, qui montre un champ de déformation plastique équivalente. Le maillage représente une petite partie (400 µm × 300 µm × 100 µm) d'une éprouvette en cuivre OFHC soumise à une traction simple, qui a subi une série de polissages *post mortem* afin de mettre en évidence la structure granulaire 3D, et qui a été maillée en respectant les joints de grains. Il y a eu pour chaque nouvelle surface une analyse EBSD, si bien qu'il a été possible d'affecter à chaque grain son orientation réelle.

Les développements les plus récents du calcul de microstructures ont pour but de caractériser de façon statistique le comportement d'un agrégat, en considérant à la fois un nombre relativement important de grains (jusqu'à plusieurs centaines), mais aussi un grand nombre d'éléments finis dans chaque grain. Il s'agit maintenant non seulement d'effectuer le passage du micro vers le macro, comme dans les premières tentatives, mais aussi d'obtenir des résultats de relocalisation afin de caractériser de façon fine les champs locaux [58–61]. Traditionnellement, les images 3D des microstructures sont générés à partir de mosaïques de Voronoï. Le plus simple est ensuite de projeter l'image sur un maillage régulier (et dans ce cas, les joints de grains traversent les éléments finis, on a des éléments « multi-phase »). Il est préférable cependant de contruire un maillage qui respecte les joints de grains, car les éléments multi-phase sont trop « raides » au niveau des joints de grains et peuvent introduire des artefacts numériques sur les champs locaux.

La Fig. 13 présente un résultat typique sur un maillage régulier (Fig. 13a), avec éléments multi-phases. On retrouve toujours pour les champs de déformation équivalente la morphologie traditionnelle qui est montrée en Fig. 13b, avec des bandes de localisation qui passent d'un grain à l'autre. Il ne s'agit pas des lignes d'instabilités que l'on retrouve avec les théories mécaniques, puisque celles-ci ne se développent pas lorsque le matériau présente un écrouissage positif. Ces bandes sont plutôt présentes pour assurer la propagation de la déformation d'un grain à l'autre, en respectant la nécessaire continuité géométrique, en l'absence d'endommagement, ce qui rend le paramètre
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Un résultat typique sur un agrégat polycristallin : (a) un maillage cube 28 × 28 × 28 comportant 200 grains ; (b) champ de déformation axiale ; (c) champ de contrainte équivalente de von Mises [62].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Courbes locales contrainte-déformation pour 200 grains, avec différentes hypothèses de calcul, pour un chargement en traction jusqu'à $e_{zz} = 0,015$. Composantes axiales : (a) modèle auto-cohérent, (b) conditions homogènes en déformation, (c) quatre faces latérales libres. Composantes latérales : (d) modèle auto-cohérent, (e) conditions homogènes en déformation, (f) quatre faces latérales libres.</figcaption></figure>

«orientation cristalline» secondaire. Au contraire, il est frappant de constater que les contours du champ de contraintes (contrainte de von Mises en Fig. 13c) respecte les joints de grains, et reste donc pleinement sensible à la morphologie cristalline.

Sans prétendre détailler tous les résultats qui sont issus de telles simulations, on donne deux informations concernant les hétérogénéités des champs locaux, car cette information doit être gardée en mémoire lorsqu'on cherche à tirer des conclusions d'observations de microstructures : le premier résultat concerne l'hétérogénéité au niveau des grains, le second l'hétérogénéité à l'intérieur d'un grain.

La Fig. 14 est construite à l'aide de données moyennes sur les grains sous un chargement de traction, avec un modèle auto-cohérent (approche de [63]) (Fig. 14a,d) et par éléments finis, avec deux sortes de conditions aux limites, traction simple (Fig. 14b,e) et déformations imposées (Fig. 14c,f). Pour chaque
<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) Déformation plastique équivalente (von Mises) *versus* distance aux joints de grains ; (b) Contrainte équivalente (von Mises) *versus* distance aux joints de grains.</figcaption></figure>

cas, on reporte les courbes de traction (Fig. 14a,b,c), et les courbes donnant la contrainte latérale en fonction de la contrainte axiale. On retrouve bien en Fig. 14a la relation simple entre déformations et contraintes locales, les contraintes résiduelles étant simplement proportionnelles à la différence entre la déformation plastique moyenne et la déformation plastique locale, ce qui se solde par une redistribution linéaire. Ce même diagramme contrainte-déformation n'a pas de construction aussi rigoureuse dans le cas des calculs par éléments finis. On note toutefois, ce qui est assez intuitif, que lorsque les déformations sont imposées aux frontières du cube, la dispersion des déformations locales dans chaque grain est relativement faible, et celle des contraintes est relativement élevée (Fig. 14b). Le contraire est vrai lorsqu'on laisse les quatre faces latérales du cube libres (Fig. 14c).

Il est remarquable de noter que les contraintes latérales sont élevées, pour les trois cas de chargement. Ceci rappelle que l'état de traction simple (une seule composante non nulle dans le tenseur de contrainte) n'existe pas dès qu'on s'intéresse à l'échelle de la microstructure. Les valeurs les plus élevées des contraintes latérales sont obtenues pour le calcul qui est réalisé à déformations imposées.

L'hétérogénéité constatée au niveau des grains est encore magnifiée lorsqu'on considère l'intérieur de ceux-ci. La Fig. 15 montre des état de contrainte (Fig. 15a) et de déformation (Fig. 15b) à l'intérieur d'un grain d'un agrégat sollicité en traction. Les valeurs sont reportées en fonction de la distance aux joints de grains. On observe que les points situés au centre du grain sont soumis aux déformations et contraintes moyennes. Au contraire, à proximité des joints, on obtient à la fois des contraintes et des déformations plus grandes et plus petites que ces valeurs moyennes.

# 6. PROBLÈMES COUPLÉS

## 6.1 Modèles de CSC

### 6.1.1 Modèle

Comme on l'a dit précédemment, les joints de grains (ou plutôt les « zones affectées par les joints de grains ») sont le siège de phénomènes particuliers à l'intérieur de l'agrégat, en raison des défauts qu'implique la présence du joint et des incompatibilités de déformation. On propose ici une approche de la corrosion sous contrainte du zircaloy (Zy4) qui est utilisé dans les centrales nucléaires de type PWR [64, 65]. Dans le process réel, les pastilles de dioxyde d'uranium contenues dans les tubes en zircaloy gonflent et peuvent sous certaines conditions combler le jeu initial, pour finalement produire une « interaction pastille-gaine » qui génère des fissures intergranulaires dont la propagation est favorisée
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Modèle géométrique d'un agrégat comportant des zones de comportement affecté par les joints de grains : (a) maillage 2D, (b) joints d'un maillage 3D, (c) zoom sur la zone affectée [64].</figcaption></figure>

par la présence d’iode. Après une propagation intergranulaire le long de quelques grains, la fissure peut devenir transgranulaire (clivage) et conduire à la rupture rapide du tube. L’effet mécanique seul ne peut pas expliquer la naissance des fissures. Il s’agit bien d’un phénomène couplé, dont la modélisation nécessite de réaliser un calcul prenant en compte mécanique et endommagement, et aussi l’avancée de l’iode dans les joints de grains. Comme l’effet mécanique est lié à un problème d’indentation du tube par un morceau de dioxyde d’uranium, la zone affectée est de petite taille, et d’importants gradients de contraintes peuvent être observés sur les premiers grains à proximité de la surface. Un modèle approprié pour expliquer quantitativement les comportements locaux doit donc prendre en compte les aspects cristallins. Le modèle numérique développé [61] introduit un calcul éléments finis couplé à couplage faible : l’adsorption de l’iode est représentée par un calcul de diffusion (où la diffusivité dépend de façon forte de l’endommagement, si bien qu’il y a un très grand gradient de concentration en avant de la pointe de fissure. Dans le calcul mécanique, la zone affectée par le joint de grain est représentée par des équations de comportement particulières. Un process catastrophique se produit dans la mesure où l’avancée de l’iode favorise l’endommagement, et où l’endommagement permet à l’iode de pénétrer plus vite dans le matériau. Des maillages 2D et 3D sont disponibles, prenant en compte le corps du grain et la zone affectée (Fig. 16).

Au niveau des joints de grains, le modèle introduit donc une zone affectée qui est partagée en deux éléments, chacun étant rattaché à un grain (Fig. 16c). Dans chaque élément « zone affectée », le comportement est pour une part hérité du grain auquel il se rattache. Si on compare avec l’élément de zone cohésive du paragraphe 3.2, la présente approche offre l’avantage d’être totalement transparente (pas de souplesse parasite du joint) tant qu’il n’y a pas d’endommagement. Par contre, il faut admettre que l’épaisseur de l’élément est une donnée qui fait partie intégrante du modèle. Trois autres variables sont présentes, l’endommagement, l’ouverture et le glissement du joint de grains. Les mécanismes d’ouverture sur un plan de vecteur normal **n** produit une pseudo-déformation portée par le produit tensoriel **n** ⊗ **n**, alors que pour le glissement du joint le tenseur est la partie symétrique de **n** ⊗ **t**, **t** étant la direction de cisaillement maximum dans le plan. L’endommagement introduit dans le modèle a pour effet d’accroître les vitesses d’ouverture et de glissement. Son évolution dépend d’une variable critique équivalente qui combine la contrainte normale et le cisaillement. Une description plus détaillée du modèle et de ses capacités est faite dans [66]. Les équations ont la forme suivante :

- Elasticité et écoulement plastique

$$ \text{décomposition de la déformation : } \dot{\varepsilon} = \dot{\varepsilon}^e + \dot{\varepsilon}^p \quad (97) $$

$$ \text{loi élastique avec endommagement isotrope : } \sigma = (1 - D)\tilde{L} : \varepsilon^e \quad (98) $$
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Simulation de la diffusion préférentielle aux joints de grains : (a) Image générale ; (b) concentration dans les zones affectées par les joints [64].</figcaption></figure>

$$ \text{ouverture et glissement : } \dot{\underline{\epsilon}}^p = \dot{\delta}\underline{n} \otimes \underline{n} + \dot{\gamma}\{\underline{n} \otimes \underline{t}\} \quad (99) $$

– Lois d’écoulement pour l’ouverture et le glissement

$$ \text{ouverture : } \dot{\delta} = \left\langle \frac{\langle \sigma_{11} \rangle / (1 - D) - R_n}{K_n} \right\rangle^{n_n} \quad (100) $$

$$ \text{glissement : } \dot{\gamma} = \left\langle \frac{|\tau|/(1 - D) - R_t}{K_t} \right\rangle^{n_t} \text{signe}(\tau) \quad (101) $$

– Evolution de l’endommagement

$$ \text{variable critique : } \sigma_D = \sqrt{\sigma^2 + \beta\tau^2} \quad (102) $$

$$ \text{\'evolution de l'endommagement : } \dot{D} = \left\langle \frac{\sigma_D - R_D}{A} \right\rangle^r (1 - D)^{-k} \quad (103) $$

### 6.1.2 Résultats

La possession de maillages «équipés» de joints de grains permet de réaliser de nombreuses investigations numériques. La Fig. 17 montre ainsi une expérience numérique de diffusion préférentielle le long des joints de grains, à partir d'une surface extérieure sur laquelle la concentration est uniforme. Le rapport entre les diffusivité au cœur du grain et dans les éléments affectés par le joint de grains est de 100 (en faveur des joints de grains).

La Fig. 18 montre quant à elle les résultats typiques qui sont obtenus sur un agrégat 2D. Les cartes (a-c) représentent la propagation intergranulaire qui est obtenue lorsqu'on résout le problème mécanique seul : les fissures s'amorcent n'importe où dans le matériau, et la propagation est un mélange de progression de fissure directe et de coalescence. Pour les planches (d-f), l'influence de l'iode est introduite du côté gauche du maillage, si bien que les joints de grains sont plus critiques à ce point. Une fissure principale se forme puis se propage. On note la présence de «branchements» à chaque point triple. Les formes finales des fissures sont tout à fait réalistes. La description complète de cette approche ([66]) montre qu'on peut retrouver les valeurs de vitesse de propagation expérimentales, et aussi que l'on peut expliquer le passage rupture intergranulaire-rupture transgranulaire.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Propagation de fissure prédite par l'évolution de l'endommagement intergranulaire : (a,b,c) cas non-couplé ; (d, e, f) cas couplé.</figcaption></figure>

## 6.2 Modèles d'oxydation

L'étude du développement des couches d'oxydes et des conséquence mécaniques de leur existence peut être menée à plusieurs échelles. Ce paragraphe donne un éclairage sur quelques approches possibles, pour l'oxydation du zircaloy et pour celle du nickel. Ces exemples sont représentatifs de deux types d'oxydation différents, anodiques ou cathodiques, qui conduisent au développement de la couche d'oxyde au niveau de la surface externe (le métal diffuse au travers de la couche d'oxyde, comme pour le système Ni-NiO), ou au niveau de l'interface entre matériau de base et oxyde (l'oxygène diffuse au travers de la couche d'oxyde). Souvent on considère que l'interface métal-oxyde est plane. Ceci simplifie considérablement le problème, car on peut alors se ramener à une géométrie 1D, et construire un modèle qui dans certains cas peut être totalement analytique. On montre cependant que le fait que l'interface soit ondulée génère des contraintes normales à la surface externe qui peuvent expliquer la desquamation des couches.

### 6.2.1 Interface plane

L'essai de flexion dissymétrique permet d'évaluer les déformations dans les couches en croissance à partir de la mesure du rayon de courbure, et, pourvu qu'on émette des hypothèses sur la loi de comportement, les contraintes. La planche 19 montre une géométrie particulière utilisée pour l'étude de l'oxydation du nickel [67]. L'apparition de l'oxyde conduit à une augmentation de volume spécifique : le « rapport de Pilling et Bedworth » (PBR) est obtenu en divisant le volume d'oxyde par le volume
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Flexion de poutre liée au développement d'une couche d'oxyde : (a) géométrie de l'éprouvette ; courbe résultante [68].</figcaption></figure>

de métal consommé pour le produire. Il est en général plus grand que 1 (1,68 pour le cuivre, 1,28 pour l'aluminium, 1,99 pour le chrome, 1,77 pour le fer, 1,52 pour le nickel, 1,56 pour le zirconium). Cette variation de volume se traduit donc par des termes sur les composantes axiales du tenseur de déformation, correspondant à une « déformation libre » qui vient s'ajouter aux différences de dilatations thermiques entre métal de base et couche d'oxyde pour créer des contraintes. Les calculs théoriques qui fournissent les valeurs du PBR ne fournissent cependant pas un renseignement important : comment cette variation de volume se distribue-t-elle sur les différentes composantes du tenseur de déformation (les deux composantes dans le plan de l'oxyde et la composante normale à l'interface) ? Les calculs qui considèrent que les trois composantes axiales sont les mêmes prévoient des niveaux de contraintes de compression dans les couches qui sont totalement irréalistes (on obtient par exemple des valeurs de l'ordre de 30 GPa avec le zirconium). On obtient par contre des valeurs réalistes en supposant que la déformation s'établit essentiellement dans la direction perpendiculaire à l'interface, et qu'il n'y a qu'une faible partie de la variation de volume qui est portée par les composantes du plan.

Dans ce cas, l'application d'une théorie élémentaire de plaques (Love-Kirchhoff) fournit les équations pour les efforts de membranes, $N_1$ et $N_2$, suivants les axes $x_1$ et $x_2$ ainsi que les moments autour de ces axes, $M_1$ et $M_2$. Les premières équations (eq.104, comportement de membrane) introduisent les déformations élastiques, différences entre déformation totale et déformation libre, cette dernière cumulant la partie thermique (température multipliée par la différence de coefficients de dilatation linéaire entre substrat et oxyde) et la partie provenant des composantes dans le plan de la matrice qui définit le changement de volume spécifique, $\varepsilon^T = \Delta\alpha T + \varepsilon^{cp}$. Les secondes (eq.105, comportement de flexion) font intervenir les rayons de courbure de la plaque, $R_1$ et $R_2$.

$$ \begin{pmatrix} N_1 \\ N_2 \end{pmatrix} = \frac{Ee}{1 - \nu^2} \begin{pmatrix} 1 & \nu \\ \nu & 1 \end{pmatrix} \begin{pmatrix} \varepsilon_1 - \varepsilon^T \\ \varepsilon_2 - \varepsilon^T \end{pmatrix} \qquad (104) $$

$$ \begin{pmatrix} M_1 \\ M_2 \end{pmatrix} = \frac{Ee^3}{12(1 - \nu^2)} \begin{pmatrix} 1 & \nu \\ \nu & 1 \end{pmatrix} \begin{pmatrix} 1/R_1 \\ 1/R_2 \end{pmatrix} \qquad (105) $$

La formule de Stoney est obtenue en supposant que $\varepsilon_1 = \varepsilon_2 \approx 0$, $R_1 = R_2 = R$, $N_1 = N_2 = N$, $M_1 = M_2 = M$, et que la couche d'oxyde n'offre aucune résistance à la flexion, si bien que les moments sont directement reliés aux efforts normaux, en introduisant l'épaisseur $e_s$ du substrat, par $M = Ne_s/2$. Elle relie le rayon de courbure à la contrainte de membrane dans l'oxyde :
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Estimation du niveau de contrainte présent dans une couche d'oxyde en présence d'un comportement viscoplastique, pour NiCr [67] : (a) stratégie d'activation progressive des éléments ; (b) relaxation à l'intérieur de la couche d'oxyde.</figcaption></figure>

$$ \frac{1}{R} = \frac{6(1 - v_s)N}{E_s e_s^2} \qquad (106) $$

En fait, une éprouvette rectangulaire ne fléchit pas avec deux rayons de courbures: l'un des deux est infini, si bien qu'on se trouve plutôt en condition de déformation dans un plan, et, en supposant que les coefficients de Poisson sont les mêmes pour le substrat et la couche déposée :

$$ \frac{1}{R} = \frac{6(1 + v)E_d e_d e^{cp}}{E_s e_s^2} \qquad (107) $$

où $e_d$ est l'épaisseur du dépôt. La Fig. 19b montre la calibration de la composante dans le plan de la déformation libre, sur une expérience de flexion d'une plaque mince en zircaloy. On obtient un bon accord entre simulation et expérience pour une déformation libre de 0,5%, soit très loin du tiers du PBR [69].

Lorsque le comportement est viscoplastique, le calcul par éléments finis fournit le moyen d'évaluer des valeurs de contraintes en régimes transitoires. La Fig. 20 montre une stratégie élémentaire, qui consiste à créer un maillage dans lequel une couche d'éléments initialement vides sont progressivement « activés » (Fig. 20a) pour rendre compte de l'avancement de l'oxydation. Cette méthode est un peu grossière, dans la mesure où elle introduit une discrétisation spatiale, donc une discontinuité de la propagation qui n'existe pas dans la réalité. Elle permet néanmoins d'avoir une estimation des contraintes qui se développent dans la couche en croissance, et leur relaxation en fonction du temps (Fig. 20b). Plus la couche d'oxyde est épaisse, moins les contraintes maximales de compression sont élevées dans ce cas là.

En résumé, l'approche macroscopique est susceptible de fournir une bonne évaluation du comportement contrainte-déformation dans la couche d'oxyde, pour les composantes parallèles à l'interface. Elle ne donne pas d'information sur les mécanismes locaux, ni sur les contraintes qui sont définies sur des facettes parallèles à la surface. Il faut donc envisager de compliquer un peu l'approche, en particulier lorsque l'interface métal-oxyde n'est pas plane.

### 6.2.2 Génération d'inter interfaces ondulées et conséquences de leur présence

La Fig. 21 montre un schéma de modélisation du champ de contraintes en présence d'une interface ondulée dans un tube de zirconium [69]. Dans la mesure où la période de l'ondulation est petite devant le
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Conséquence de la géométrie de l'interface sur la modélisation du champ de contraintes locales dans le cas du système Zr/ZrO₂ : (a) simulation de l'ondulation à l'aide d'un maillage en conditions axisymétriques ; (b) distribution des contraintes de compression et de traction.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Champs locaux dans une couche d'oxyde pour une interface ondulée : (a) contrainte dans le plan ; (b) contrainte normale au plan ; (c) contrainte équivalente de von Mises ; (d) déformation plastique équivalente de von Mises [68].</figcaption></figure>

rayon du tube, on peut l'isoler dans un petit cylindre dont le rayon fait quelques dizaines de micromètres, et effectuer le calcul en axisymétrique, l'axe du cylindre étant porté par le rayon du tube (Fig. 21a). La zone où la couche d'oxyde est la plus épaisse a subi un important gonflement radial, lié au PBR, puisque la déformation libre perpendiculairement à l'interface accommode pratiquement toute la variation de volume, soit plus de 50%. Bridée dans ce gonflement par la zone qui l'entoure (où l'oxyde est mince), elle est donc en compression, alors que la zone mince est en traction. Cet état de contrainte est le résultat d'un auto-équilibre entre les deux zones. Un point expérimental vient confirmer l'intérêt de cette analyse : les fissures qui se développent dans la couche d'oxyde se situent justement aux endroits où la couche d'oxyde est la plus mince. Si au contraire on s'intéresse à la contrainte dans le plan, cet auto-équilibre a lieu entre la couche d'oxyde et le substrat : ceci conduit donc à une compression dans la couche d'oxyde et une traction dans le substrat (qui reste faible, en raison de la grande épaisseur du substrat devant celle de la couche d'oxyde) (Fig. 21b).

Les cartes des différents champs locaux sont reportés en Fig. 22. On note que la couche d'oxyde est le siège des contraintes les plus importantes, et que le niveau de contrainte induit par la couche dans le substrat suffit à déclencher localement de l'écoulement plastique [68].

Le développement de l'ondulation de l'interface peut être expliqué de plusieurs manières. Un calcul sur une géométrie plane a permis de formuler une hypothèse, basée sur l'hétérogénéité du coefficient
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Distribution aléatoire du coefficient de diffusion pour un calcul plan : (a) variation en fonction de la concentration d'oxygène ; (b) champ aléatoire 2D ; (c) champ aléatoire 1D.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 24.** Contraintes résultant de l'oxydation progressive du zircaloy : (a) contrainte normale dans le plan de l'oxyde ; (b) contrainte perpendiculaire à la surface.</figcaption></figure>

de diffusion [70]. Le système étudié est un carré, initialement formé de métal de base uniquement. On applique une concentration d'oxygène donnée sur l'arête inférieure, et on réalise un calcul de diffusion. On suppose que la transformation métal-oxyde s'effectue lorsque la concentration d'oxygène atteint la composition stoechiométrique (0,59). La transition se manifeste par le fait que le coefficient de diffusion dépend de la concentration d'oxygène. La Fig. 23a montre quelques courbes typiques qui ont été retenues, et qui retiennent une valeur relativement faible dans le métal de base (concentration d'oxygène inférieure à 0,5), et qui augmente de façon aléatoire lorsqu'on atteint les concentrations présentes dans l'oxyde. Un premier calcul, réalisé avec un coefficient de diffusion totalement aléatoire dans l'espace (Fig. 23b), prévoit que le front reste plan. Au contraire, si l'aspect aléatoire n'est présent que dans une seule direction, (Fig. 23c), les propriétés sont uniforme le long de « tubes » perpendiculaires à la surface extérieure. On peut voir dans ce modèle une représentation des grains colonnaires qui constituent l'oxyde. Le front présente alors une ondulation, liée à une auto-organisation locale avec une portée d'une dizaine de tubes. Les états de contraintes qui en résultent sont reportés en Fig.24. On retrouve bien les
mêmes ordres de grandeur que pour le calcul axisymétrique. La contrainte dans le plan est uniformément en compression dans la couche d'oxyde. La contrainte normale à la surface est bien entendu nulle sur la surface libre (arête inférieure du carré), mais elle se développe de façon très locale au voisinage de l'interface, et atteint des niveaux importants en traction dans les zones où la couche d'oxyde est la plus fine. Ces niveaux de contrainte suffisent à expliquer l'apparition de fissures dans la couche d'oxyde.

# 7. CONCLUSION ET PERSPECTIVES

La plasticité cristalline est maintenant une théorie suffisamment établie pour être utilisée de façon opérationnelle dans un grand nombre d'applications industrielles (sur composants monocristallins en particulier) ou plus académiques (le calcul de microstructures pour caractériser le comportement global et local des matériaux hétérogènes). Le but de cette partie était de passer rapidement en revue les possibilités des modèles, en mentionnant leurs aspects pratiques et leurs liens avec les approches classiques de la plasticité et de la viscoplasticité. Une série de sujets n'ont pas été abordés ici, mais qui seraient utiles à la poursuite du dialogue matériau-mécanique.

L'une des extensions principales du cadre classique est proposée par la théorie des milieux continus généralisés. Le lecteur concerné pourra lire une synthèse récente dans [71]. Ces théories ont déjà été appliquées au monocristal (voir par exemple [72, 73]). Ce type d'approche permet en particulier de prendre en compte de façon explicite la taille de l'élément de volume. Il faut aussi considérer les connexions avec d'autres types d'approches, notamment les approches capables de descendre à des échelles submicroniques, et de reconnaître le caractère discret des mécanismes de déformation. C'est le cas de la dynamique discrète des dislocations (DDD), et de la dynamique moléculaire (DM), qui ont chacune leurs propres échelles en temps et en espace. Dans un futur proche, leurs progrès vont permettre de nourrir de façon pertinente les modèles des échelles plus grossières, comme ceux de plasticité cristalline : on aura à disposition des chaînes de calcul traversant les échelles, où la DDD prendra ses informations de la DM, et la plasticité cristalline de la DDD.

Si on ne se restreint pas à la mécanique, et qu'on veut aborder – comme on l'a fait en fin de chapitre – les problèmes multiphysiques, on doit compter sur les progrès des méthodes de champs de phases, qui ne seront plus confinées dans des codes ultra-spécialisées, mais vont devenir disponibles dans des codes de calcul généralistes, ce qui leur ouvrira des possibilités d'emploi nouvelles, en couplage avec d'autres modèles, notamment la mécanique. Le lecteur pourra consulter [74] dans ce recueil pour une présentation d'une implémentation dans un code de calcul par éléments finis.

Pour conclure, il faut insister sur le fait que, dans le domaine de la modélisation des matériaux plus qu'ailleurs, la modélisation numérique ne se développe pas aux dépens des aspects expérimentaux, mais qu'au contraire elle se nourrit de ces résultats. Les progrès récents, et ceux à venir, qui touchent les mesures de champs, en 2D et en 3D, offrent une occasion unique d'effectuer des comparaisons de plus en plus documentées entre modèles numériques et expériences. On peut par exemple bénéficier de mesure de champs précises en surface de spécimens [57, 75], mais également d'informations 3D sur les microstructures et les populations de fissures [76].

# Remerciements

Le document contient des planches issues de plusieurs thèses récentes, notamment celles de Fabienne Hanriot et de Sylvain Flouriot, soutenues par Snecma Moteurs ; celle de Marion Parise, soutenue par le CEA Grenoble ; celles de Olivier Diard puis Andrey Musienko, soutenues par EDF et le CEA. De nombreux apports, directs ou indirects, proviennent également du travail réalisé par Fabrice Barbe, avec la collaboration de Samuel Forest. Que tous soient remerciés pour leurs discussions vivifiantes.
Références

[1] Halphen B. and Nguyen Q.S. *J. de Mécanique*, **14** (1975) 39.

[2] Germain P., Nguyen Q.S. and Suquet P. *J. Appl. Mech.*, **50** (1983) 1010.

[3] Lemaitre J. and Chaboche J.-L. Mechanics of Solid Materials. (Cambridge University Press, Cambridge, U.K., 1990).

[4] Besson J., Le Riche R., Foerch R. and Cailletaud G. *Revue Européenne des Éléments Finis*, **7** (1998) 567.

[5] Koiter W.T. *General theorems for elastic-plastic solids*, volume 6, pages 167–221. North-Holland Publishing Company, 1960.

[6] Mandel J. *Int. J. Sol. Structures*, **1** (1965) 273.

[7] Méric L., Poubanne P. and Cailletaud G. *J. of Engng. Mat. Technol.*, **113** (1991) 162.

[8] Hanriot F., Cailletaud G. and Rémy L. In A.D. Freed and K.P. Walker, editors, *High Temperature Constitutive Modeling – Theory and Application*. ASME, 1991.

[9] Taylor G.I. *J. Inst. Metals.*, **62** (1938) 307.

[10] Bishop J.F.W. and Hill R. *Philosophical Magazine*, **42** (1951) 414.

[11] Chin G.Y. and Mamme W.L. *Trans. of the Met. Soc. of AIME*, **245** (1969) 1211.

[12] Jordan E.H. and Walker K.P. In K. Miller, editor, *Conf. on Multiaxial Fatigue*, Sheffield, (1985).

[13] Busso E.P. and McClintock F.A. *Int. J. of Plasticity*, **12** (1996) 1.

[14] Kocks U.F. and Brown T.J. *Acta Metall.*, **14** (1966) 87.

[15] Kocks U.F. *Metallurgical Transaction*, **1** (1970) 1121.

[16] Franciosi P. *Acta Metall.*, **33** (1985) 1601.

[17] Franciosi P. *Revue de Physique Appliquée*, **33** (1985) 1601.

[18] Devincre B., Kubin L. and Hoc T. *Scripta Metallurgica and Materiala*, **54** (2006) 741.

[19] Asaro R.J. *J. of Applied Mechanics*, **50** (1983) 921.

[20] Tabourot L., Fivel M. and Rauch E. *Material Science and Engineering*, A234-236 (1997) 639.

[21] Pierce D., Asaro R.J. and Needleman A. *Acta Metall.*, **31** (1985) 1951.

[22] Asaro R.J. and Needleman A. *Acta Metall.*, **33** (1985) 923.

[23] Méric L. and Cailletaud G. *J. of Engng. Mat. Technol.*, **113** (1991) 171.

[24] Cuitino A.M. and Ortiz M. *Modelling Simul. Mater. Sci. Eng.*, **1** (1992) 225.

[25] Anand L. and Kothari M.J. *Mech. Phys. Sol.*, **44** (1996) 525.

[26] Schröder J. and Miehe C. *Computational Materials Science*, **9** (1997) 168.

[27] Simo J.C. and Hughes T.R.J. *Computational Inelasticity*. Springer Verlag, 1997.

[28] McGinty R.D. and McDowell D.L. *Int. J. of Plasticity*, **22** (2006) 996.

[29] Foerch R., Besson J., Cailletaud G. and Pilvin P. *Comp. Meth. Appl. Mech. Engng*, **141** (1997) 355.

[30] Dugdale D.S. *J. Mech. Phys. Sol.*, **8** (1960) 100.

[31] Barenblatt G.I. *Advances in Appl. Mech.*, **7** (1962) 55.

[32] Needleman A. *J. of Applied Mechanics*, **54** (1987) 525.

[33] Suquet P. In J.-J. Moreau and P.D. Panagiotopoulos, editors, *Nonsmooth mechanics and applications*, pages 280–340, Springer, 1987.

[34] Michel J.-C., Suquet P. and Thébaud F. *Revue Européenne des Éléments Finis*, **3** (1994) 573.

[35] Lebon F., Rizzoni R. and Ronel-Idrissi S. *Computers and Structures*, **82** (2004) 1929.

[36] Tvergaard V. *Material Science and Engineering*, A125 (1990) 203.

[37] Allix O. and Ladevèze P. *J. Compos. Struct.*, **22** (1992) 235.

[38] Chaboche J.-L., Girard R. and Schaff A. *Computational Mechanics*, **20** (1997) 3.

[39] Chaboche J.-L., Feyel F. and Monerie Y. *Int. J. Sol. Structures*, **38** (2001) 3127.

[40] Gallerneau F. *Etude et modélisation de l'endommagement d'un superalliage monocristallin revêtu pour aube de turbine*. PhD thesis, École Nationale Supérieure des Mines de Paris, 1995.
[41] Cailletaud G., Chaboche J.-L., Forest S. and Rémy L. *Revue de Métallurgie*, février 2003 (2001) 165.

[42] Nouailhas D., Culié J.-P., Cailletaud G. and Méric L. *Eur. J. Mech., A/Solids*, **14** (1994) 137.

[43] Flouriot S., Forest S., Cailletaud G., Köster A., Rémy L., Burgardt B., Gros V. and Mosset J., Delautre S. *Int. J. Frac*, **124** (2003) 43.

[44] Rice J. *Mech. of Materials*, **6** (1987) 317.

[45] Sai K., Cailletaud G. and Forest S. *Mech. of Materials*, **38** (2006) 203.

[46] Engel J.-J. *Modélisation du comportement global des métaux en plasticité et en viscoplasticité*. PhD thesis, École Nationale Supérieure des Mines de Paris, 1978.

[47] Kalidindi S.R., Bronkhorst C.A. and Anand L. *J. Mech. Phys. Sol.*, **40** (1992) 536.

[48] Miehe C., Schröder J. and Schotte J. *Comp. Meth. Appl. Mech. Engng*, **171** (1999) 387.

[49] Bugat S., Besson J. and Pineau A. *Computational Materials Science*, **16** (1999) 158.

[50] Mika D.P. and Dawson P.R. *Acta Metall.*, **47** (1999) 1355.

[51] Evers L.P., Parks D.M., Brekelmans W.A.M. and Geers M.G.D. *J. Mech. Phys. Sol.*, **50** (2002) 2403.

[52] Kim D.J., Ku T.W. and Kang B.S. *J. of Materials Processing Technology*, **130** (2002) 456.

[53] Delaire F., Raphanel J.L. and Rey C. *Acta Metall.*, **48** (2000) 1075.

[54] Raabe D., Sachtleber M., Zhao Z., Roters F. and Zaefferer S. *Acta Metall.*, **49** (1981) 3433.

[55] Zeghadi A., N’guyen F., Forest S., Gourgues A.-F. and Bouaziz O. *Philosophical Magazine*, **87** (2007) 1401.

[56] Zeghadi A., Forest S., Gourgues A.-F. and Bouaziz O. *Philosophical Magazine*, **87** (2007) 1425.

[57] Musienko A., Tatschl A., Schmidegg K., Kolednik O., Pippan R. and Cailletaud G. *Acta Mat.*, **55** (2007) 4121.

[58] Barbe F., Decker L., Jeulin D. and Cailletaud G. *Int. J. of Plasticity*, **17** (2001) 513.

[59] Barbe F., Forest S. and Cailletaud G. *Int. J. of Plasticity*, **17** (2001) 537.

[60] Diard O., Leclercq S., Rousselier G. and Cailletaud G. *Computational Materials Science*, **25** (2002) 73.

[61] Cailletaud G., Diard O. and Musienko A. Ahzi In S., Khaleel M.A., Zbib H.M., Zikry M.A., LaMatina B. and Cherkaoui M., editors, *Multiscale Modelling and Characterisation of Elastic-Inelastic Behaviour of Engineering Materials (Proc.of the IUTAM Symposium on Multiscale Modelling of Engng. Materials. Marrakech, Oct. 2002)*, pages 149–156. Kluwer, 2004.

[62] Diard O., Leclercq S., Rousselier G. and Cailletaud G. *Int. J. of Plasticity*, **21** (2005) 691.

[63] Berveiller M. and Zaoui A. *J. Mech. Phys. Sol.*, **26** (1979) 325.

[64] Musienko A. *Plasticité cristalline en présence de grandes déformations et d'endommagement*. PhD thesis, École Nationale Supérieure des Mines de Paris, 19 mars 2005.

[65] Diard O. *Un exemple de couplage comportement-endommagement-environnement dans les polycristaux. Application à l'interaction pastille-gaine*. PhD thesis, École Nationale Supérieure des Mines de Paris, 2001.

[66] Musienko A. and Cailletaud G. *Soumis*, 2009.

[67] Huntz A.-M., Calvarin-Amiri G., Evans H.E. and Cailletaud G. *Oxidation of Metals*, **57** (2002) 499.

[68] Parise M., Sicardy O. and Cailletaud G. *J. Nucl. Mat.*, **256** (1998) 35.

[69] Parise M. *Mécanismes de corrosion des alliages de zirconium : étude des cinétiques initiales d'oxydation et du comportement du système métal-oxyde*. PhD thesis, École Nationale Supérieure des Mines de Paris, 1996.

[70] Parise M., Foerch R. and Cailletaud G. *Journal de Physique IV*, **256** (1999) 311.

[71] Papenfuss C. and Forest S. *J. of Non-Equilibrium Thermodynamics*, **31** (2006) 319.

[72] Forest S., Barbe F. and Cailletaud G. *Int. J. Sol. Structures*, **37** (2000) 7105.

[73] Forest S., Sievert R. and Aifantis E.C. *J. of the Mechanical Behaviour of Materials*, **13** (2002) 219.
[74] Ammar K., Appolaire B., Cailletaud G., Feyel F. and Forest S. In *Mécanismes et Mécanique des Interactions Plasticité - Environnement : Oxydation, Corrosion, Hydrogène*, Argelès-sur-Mer, 2009. EDP Sciences.

[75] Kempf D., Vignal V., Cailletaud G., Oltra R., Weeber J.C. and Finot E. *Philosophical Magazine*, **87** (2007) 1379.

[76] Buffière J., Ferrié E., Proudhon H. and Ludwig W. *Material Science and Technology*, **22** (2006) 1019.