# Aspects numériques d'un modèle variationnel de plasticité cristalline

Laurent Stainier

► To cite this version:

Laurent Stainier. Aspects numériques d'un modèle variationnel de plasticité cristalline. 9e Colloque national en calcul des structures, CSMA; LMT, May 2009, Giens, France. hal-01400483

HAL Id: hal-01400483

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

Submitted on 22 Nov 2016

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
# Aspects numériques d'un modèle variationnel de plasticité cristalline

L. Stainier

Institut de recherche en Génie Civil et Mécanique (GeM - UMR 6183 CNRS)
Ecole Centrale de Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes cedex 3
Laurent.Stainier@ec-nantes.fr

**Résumé** — Cet article décrit l'application d'une formulation variationnelle de la mécanique au cas de modèles de comportement de type plasticité cristalline, formulation offrant une solution intéressante au problème de la détermination de l'ensemble des systèmes de glissement actifs. Nous explorons différentes options possibles pour intégrer des modèles existant dans ce cadre, en mettant en évidence les conséquences de ces choix, à la fois en terme de précision et de robustesse de l'algorithme numérique. Les aspects thermo-mécaniques sont également brièvement abordés.

**Mots clés** — plasticité cristalline, formulation variationnelle, modèles constitutifs.

## 1 Introduction

Le comportement mécanique macroscopique des métaux polycristallins est régi par les micro-mécanismes à l'oeuvre à l'échelle des grains individuels et de leurs interfaces. Parallèlement aux moyens modernes d'observation expérimentale (microscopie électronique), la simulation numérique directe du comportement d'agrégats polycristallins offre une possibilité complémentaire de comprendre ces micromécanismes, au moins de manière qualitative. On peut également espérer que ce type d'approche puisse conduire à une modélisation prédictive, par homogénéisation numérique, du comportement des matériaux polycristallins, même si cet objectif n'est encore que très partialement atteint.

Le comportement des grains monocristallins peut être décrit par des modèles de plasticité cristalline, représentant l'effet statistique du mouvement des dislocations. La mise en oeuvre algorithmique de ces modèles de plasticité cristalline dans le cadre incrémental propre aux simulations numériques non-linéaires pose cependant encore un certain nombre de problèmes : la détermination des systèmes de glissement actifs, spécialement dans le cas de chemins de chargement non proportionnels, peut poser problème. En effet, la complexité des modèles (on peut compter jusqu'à 48 systèmes de glissement potentiels), combinée à la nature incrémentale de l'approche, peut conduire à des problèmes d'unicité de la solution numérique. Les algorithmes utilisés font alors partie intégrante du modèle, et le choix d'un algorithme de résolution correspond donc au choix d'un critère permettant de déterminer la "bonne" solution parmi les solutions au problème mathématique. La formulation variationnelle que nous avons proposée en [5] offre de ce point de vue une option attrayante : la solution choisie est celle qui minimise une énergie incrémentale directement reliée au travail extérieur. Pour comparaison, l'algorithme proposé par Anand [1] choisit la solution qui minimise l'incrément total de glissement. On peut également mentionner l'algorithme proposé par Cuitiño [3], qui présente beaucoup de points communs avec le nôtre. On trouvera également une discussion d'autres algorithmes dans [2]. A ces algorithmes implicites, on
peut ajouter les algorithmes explicites [4]. Ceux-ci considèrent comme actifs tous les systèmes de glissement pour lesquel la valeur critique de la cission (contrainte projetée sur le système de glissement) est atteinte au début du pas de temps, et déterminent l'amplitude des glissements plastiques à partir d'une relation cinétique. Le recours à cette relation limite cependant l'applicabilité des algorithmes explicites aux modèles constitutifs dépendant de la vitesse de déformation (visco-plasticité), ce qui n'est pas le cas des algorithmes implicites. De plus, les algorithmes explicites ne donnent des résultats satisfaisant que pour des pas de temps relativement petits.

# 2 Formulation variationnelle

## 2.1 Équations du problème

Nous nous plaçons dans le cadre d'une approche locale de la mécanique des solides, où l'état thermo-mécanique en un point matériel est entièrement représenté par le gradient de déformation $F$, la température absolue $T$, et un ensemble de variables internes à définir. Les solides visco-plastiques sont caractérisés par l'existence d'une classe de déformations $F^p$, dites déformations plastiques, caractérisées par une absence de distorsion ou de rotation du réseau cristallin, et ne causant par conséquence aucune contrainte à longue portée. En plus des déformations plastiques, on peut supposer l'existence d'un certain niveau de distorsion du réseau cristallin, correspondant à des déformations réversibles $F^e$, dites élastiques. On peut alors écrire :

$$F = F^e F^p$$

Les déformations plastiques dans les monocristaux sont cristallographiques par nature, et, suivant Rice [7], on peut écrire la loi d'écoulement plastique suivante :

$$\dot{F}^p = \left( \sum_{\alpha=1}^{N} \dot{\gamma}^{(\alpha)} s^{(\alpha)} \otimes m^{(\alpha)} \right) F^p \quad (1)$$

où $\dot{\gamma}^{(\alpha)}$ est le glissement plastique, et $s^{(\alpha)}$ et $m^{(\alpha)}$ sont les vecteurs unitaires définissant la direction et le plan de glissement relatifs au système de glissement $\alpha$. On considérera ici que la collection $\gamma$ des glissements plastiques constitue avec la déformation plastique $F^p$ l'ensemble des variables internes du modèle.

On postule alors l'existence d'une énergie libre de Helmholtz de la forme suivante :

$$W(F, F^p, \gamma, T) = W^e(C^e, T) + W^p(\gamma, T)$$

où $C^e = F^e T F^e$ est le tenseur élastique de Cauchy-Green à droite. Le tenseur des contraintes de Piola $P$ dérive de l'énergie libre

$$P = \frac{\partial W}{\partial F} = F^e \left( 2 \frac{\partial W^e}{\partial C^e} \right) F^{p-T}$$

Considérant la loi d'écoulement plastique (1) comme une liaison non-holonome entre la déformation plastique $F^p$ et les glissements plastiques $\gamma$, la force thermodynamique conjuguée à ceux-ci est donnée par

$$Y^{(\alpha)} = -\frac{\partial W}{\partial \gamma^{(\alpha)}} = (F^e T F^{pT}) \cdot (s^{(\alpha)} \otimes m^{(\alpha)}) - \frac{\partial W^p}{\partial \gamma^{(\alpha)}} = \tau \cdot [(F^e s^{(\alpha)}) \otimes (m^{(\alpha)} F^{e-1})] - g^{(\alpha)} \quad (2)$$

où $\tau = PF^T$ est le tenseur des contraintes de Kirchhoff, et $g^{(\alpha)} = \partial W^p / \partial \gamma^{(\alpha)}$ une contrainte limite d'écoulement plastique correspondant à des mécanismes de stockage (irréversible) d'énergie. On voit également apparaître dans (2) la cission $\tau^{(\alpha)}$ :

$$\tau^{(\alpha)} = \tau \cdot [(F^e s^{(\alpha)}) \otimes (m^{(\alpha)} F^{e-1})]$$
De façon très générale, on peut formuler les lois cinétiques d'évolution des glissements plastiques au travers d'un potentiel convexe de dissipation $\psi(Y)$ :

$$ \dot{\gamma}^{(\alpha)} = \frac{\partial \psi(Y)}{\partial Y^{(\alpha)}} $$

ou de son dual (convexe lui aussi)

$$ \psi^*(\dot{\gamma}) = \sup_Y [Y \cdot \dot{\gamma} - \psi(Y)] \quad ; \quad Y^{(\alpha)} = \frac{\partial \psi^*(\dot{\gamma})}{\partial \dot{\gamma}^{(\alpha)}} $$

Les équations locales régissant l'écoulement plastique s'écrivent donc

$$ \tau^{(\alpha)} = g^{(\alpha)}(\gamma, T) + \psi_{,\dot{\gamma}^{(\alpha)}}^*(\dot{\gamma}; \gamma, T) \quad (3) $$

où on a explicité le fait que le potentiel de dissipation pouvait dépendre paramétriquement de l'état courant du matériau.

## 2.2 Formulation incrémentale

En vue de l'utilisation de modèles de plasticité cristalline dans le cadre de calculs de structures, un problème fondamental est celui de l'intégration numérique des équations (3) sur un incrément de temps (et/ou de chargement) fini. L'approche variationnelle proposée en [5] présente de ce point de vue une solution intéressante.

Considérons un pas de temps générique $[t_n, t_{n+1}]$, pour lequel l'état initial $\{F_n, F_n^p, \gamma_n\}$ et la déformation $F_{n+1}$ sont connus. La première étape pour passer des équations continues aux équations incrémentales est d'établir une règle de mise à jour de $F^p$ compatible avec la loi d'écoulement continue (1) :

$$ F_{n+1}^p = \exp \left[ \sum_{\alpha=1}^{N} \Delta\gamma^{(\alpha)} s^{(\alpha)} \otimes m^{(\alpha)} \right] F_n^p $$

avec $\Delta\gamma = \gamma_{n+1} - \gamma_n$. L'utilisation d'une exponentielle tensorielle dans le cadre d'une loi d'écoulement plastique incrémentale a initialement été proposée dans [10].

Les glissements plastiques en fin de pas sont alors obtenus via le problème de minimisation suivant :

$$ \mathcal{W}^{\text{eff}}(F_{n+1}; F_n, F_n^p, \gamma_n) = \inf_{\gamma_{n+1}} \mathcal{W}(F_{n+1}, F_{n+1}^p, \gamma_{n+1}; F_n, F_n^p, \gamma_n) \quad (4) $$

où la densité d'énergie incrémentale $\mathcal{W}$ est donnée par

$$ \mathcal{W}(F_{n+1}, F_{n+1}^p, \gamma_{n+1}; F_n, F_n^p, \gamma_n) = W(F_{n+1}, F_{n+1}^p, \gamma_{n+1}) - W(F_n, F_n^p, \gamma_n) + \Delta t \Psi^*\left(\frac{\Delta\gamma}{\Delta t}; \gamma_{n+\epsilon}\right) $$

avec $\gamma_{n+\epsilon} = (1 - \epsilon)\gamma_n + \epsilon\gamma_{n+1}$ ($\epsilon \in [0, 1]$). On observe en effet que la condition de stationnarité associée au problème de minimum (4) s'écrit

$$ \tau_{n+1}^{(\alpha)} = g^{(\alpha)}(\gamma_{n+1}) + \psi_{,\dot{\gamma}^{(\alpha)}}^*(\frac{\Delta\gamma}{\Delta t}; \gamma_{n+\epsilon}) + \epsilon\Delta t \Psi_{,\dot{\gamma}^{(\alpha)}}^*(\frac{\Delta\gamma}{\Delta t}; \gamma_{n+\epsilon}) \quad (5) $$

On peut noter que la dépendance à la température a été omise dans la relation ci-dessus pour des raisons de clarté. On constate bien que la relation incrémentale (5) est consistante avec la relation continue (3) pour $\Delta t \to 0$. En pratique, la solution du problème de minimum (4) doit également vérifier la contrainte d'inégalité $\gamma_{n+1}^{(\alpha)} \ge \gamma_n^{(\alpha)}$ ($\forall \alpha = 1, \dots, N$). Cette contrainte peut éventuellement être implicitement incluse via une expression non régulière de $\psi^*$, mais cette option n'est pas nécessairement pratique d'un point de vue algorithmique, et on préférera traiter le problème comme une minimisation sous contraintes (de borne).

Comme montré en [5], la densité d'énergie incrémentale $\mathcal{W}^{\text{eff}}(F_{n+1})$ résultant de la minimisation (4) constitue un potentiel pseudo-hyperélastique pour les contraintes. Le problème aux conditions limites global peut alors être écrit à chaque pas de temps sous forme variationnelle (de type minimisation de l'énergie potentielle).
# 3 Plasticité cristalline

## 3.1 Modèles d'écrouissage

Les modèles micromécaniques d'écrouissage utilisés en plasticité cristalline sont généralement écrits moyennant l'expression d'une matrice d'écrouissage $h^{\alpha\beta}$, telle que :

$$ \dot{g}^{(\alpha)} = \sum_{\beta=1}^{N} h^{\alpha\beta} \dot{\gamma}^{(\beta)} \qquad (6) $$

Cette forme de la relation d'écrouissage résulte de l'observation de ce que la limite d'écoulement plastique $g^{(\alpha)}$ sur un système de glissement est influencée à la fois par l'activité des dislocations dans le plan de glissement en question, mais aussi par l'activité dans les autres systèmes de glissement (création d'obstacles de type forêt de dislocations). Ce second type d'écrouissage est appelé écrouissage latent.

Les modèles les plus simples (p.e. [6]) sont du type suivant :

$$ h^{\alpha\beta} = h(\gamma)(q + (1 - q)\delta_{\alpha\beta}) $$

avec $q > 1$ un coefficient constant d'écrouissage latent. Des modèles basés sur des considérations micromécaniques (analyse statistique du mouvement des dislocations) ont par ailleurs été proposés ([3], [9]), conduisant à une matrice d'écrouissage diagonale :

$$ h^{\alpha\beta} = h^{\alpha}(\gamma, g^{(\alpha)}) \delta_{\alpha\beta} $$

Malgré sa forme diagonale, ce type de modèle prend en compte l'écrouissage latent au travers de densités d'obstacles variables, elles-mêmes calculées à partir des glissements plastiques.

## 3.2 Implémentation

Repartant de la définition de la contrainte d'écoulement plastique correspondant à des mécanismes de stockage d'énergie, on peut écrire

$$ dg^{(\alpha)} = \sum_{\beta=1}^{N} \frac{\partial^2 W^p(\gamma)}{\partial\gamma^{(\alpha)} \partial\gamma^{(\beta)}} d\gamma^{(\beta)} \qquad (7) $$

Comparant cette relation (7) à (6), on est amené à identifier la matrice hessienne associée à $W^p$ à une matrice d'écrouissage. Malheureusement, les matrices d'écrouissage proposées dans la littérature sont rarement intégrables, et il n'est donc pas possible de directement leur associer une densité d'énergie plastique.

Une première solution, initialement proposée dans [5], consiste à utiliser une approximation quadratique incrémentale de la densité d'énergie plastique :

$$ W^p(\gamma_{n+1}) \approx W^p(\gamma_n) + \sum_{\alpha=1}^{N} g_n^{(\alpha)} \Delta\gamma^{(\alpha)} + \frac{1}{2} \sum_{\alpha=1}^{N} h^{\alpha}(\gamma_n, g_n^{(\alpha)}) (\Delta\gamma^{(\alpha)})^2 $$

On obtient alors une approximation explicite de la contrainte d'écoulement plastique en fin de pas :

$$ g_{n+1}^{(\alpha)} = g^{(\alpha)}(\gamma_{n+1}) = g_n^{(\alpha)} + h^{\alpha}(\gamma_n, g_n^{(\alpha)}) \Delta\gamma^{(\alpha)} \qquad (8) $$

Cette approximation permet l'utilisation d'algorithmes de minimisation sous contraintes pour déterminer l'ensemble des systèmes de glissement actifs (ceux pour lesquels $\Delta\gamma^{(\alpha)} > 0$). Cependant les modèles d'écrouissage proposés en [3] et [9] présentent un coefficient d'écrouissage potentiellement très élevé au départ (pour les petits glissements plastiques) et chutant assez rapidement
avec la déformation plastique. L'approximation explicite (8) oblige donc à faire de très petits pas en début de calcul, sous peine de surestimer significativement la contrainte d'écoulement plastique sur les systèmes actifs.

Dans le cas d'un comportement viscoplastique, la limite d'écoulement plastique se décompose additivement en une partie indépendante de la vitesse de déformation et une partie dépendante de celle-ci. Dans la première approche ci-dessus, on considère que la partie indépendante de la vitesse de déformation est entièrement associée à des mécanismes de stockage irréversibles (et donc à $W^p$), tandis que le potentiel de dissipation $\psi^*$ correspond à la partie dépendant de la vitesse de déformation. Le potentiel de dissipation prend alors, par exemple, la forme suivante :

$$ \psi^*(\dot{\gamma}; \gamma) = \begin{cases} \sum_{\alpha=1}^{N} \frac{m\tau_0\dot{\gamma}_0}{m+1} \left(\frac{\dot{\gamma}^{(\alpha)}}{\dot{\gamma}_0}\right)^{\frac{m+1}{m}} & \text{si } \dot{\gamma}^{(\alpha)} \ge 0 \quad \forall \alpha \\ +\infty & \text{si } \dot{\gamma}^{(\alpha)} < 0 \text{ pour un système au moins} \end{cases} $$

On retrouve un comportement indépendant de la vitesse de chargement pour $m \to +\infty$. On peut cependant remarquer que, d'un point de vue purement mécanique, seule compte la somme $W_{\gamma^{(\alpha)}}^p + \psi_{\gamma^{(\alpha)}}^*$, puisque c'est sous cette forme que les potentiels interviennent dans les équations de stationnarité (3). On peut donc alternativement inclure la partie écrouissable de la limite d'écoulement plastique dans le potentiel de dissipation, qui prend alors la forme

$$ \psi^*(\dot{\gamma}; \gamma) = \sum_{\alpha=1}^{N} \tau_v^{(\alpha)}(\gamma) \dot{\gamma}^{(\alpha)} + \sum_{\alpha=1}^{N} \frac{m\tau_0\dot{\gamma}_0}{m+1} \left(\frac{\dot{\gamma}^{(\alpha)}}{\dot{\gamma}_0}\right)^{\frac{m+1}{m}} $$

où la cission critique $\tau_v^{(\alpha)}$ est définie sous forme différentielle :

$$ \frac{\partial \tau_v^{(\alpha)}}{\partial \gamma^{(\alpha)}} = h^{\alpha}(\gamma, \tau_v^{(\alpha)}) \quad (\tau_v^{(\alpha)}(0) = \tau_{v0}^{(\alpha)}) \qquad (9) $$

La densité d'énergie plastique peut alors être choisie identiquement nulle (par exemple).

La forme différentielle (9) se prête très bien à une intégration numérique implicite :

$$ \tau_v^{(\alpha)}(\gamma_{n+1}) \equiv g_{n+1}^{(\alpha)} : g_{n+1}^{(\alpha)} = g_n^{(\alpha)} + h^{\alpha}(\gamma_{n+1}, g_{n+1}^{(\alpha)}) \Delta\gamma^{(\alpha)} $$

On note que, au vu de la relation incrémentale (5), il est également nécessaire de calculer

$$ \frac{d\tau_v^{(\alpha)}}{d\gamma^{(\beta)}} = \frac{h^{\alpha}(\gamma, \tau_v^{(\alpha)}) \delta_{\alpha\beta} + \frac{\partial h^{\alpha}}{\partial \gamma^{(\beta)}}(\gamma, \tau_v^{(\alpha)}) \Delta\gamma^{(\alpha)}}{1 - \frac{\partial h^{\alpha}}{\partial g^{(\alpha)}}(\gamma, \tau_v^{(\alpha)}) \Delta\gamma^{(\alpha)}} $$

Cette seconde approche implicite permet d'utiliser de plus grands incréments de chargement en début de calcul, tout en s'assurant d'une bonne approximation de l'écrouissage initial, surtout dans les cas présentant un fort écrouissage latent. Cela conduit également à une prédiction plus robuste de l'ensemble des systèmes de glissement actifs.

## 4 Aspects thermomécaniques

Si, d'un point de vue purement mécanique, la décomposition de la limite d'écoulement plastique en une partie associée à des mécanismes de stockage et une partie associée à des mécanismes dissipatifs est relativement arbitraire, il n'en est évidemment pas de même d'un point de vue thermo-mécanique. En effet, la quantité de chaleur dégagée lors d'un processus de déformation plastique est directement liée à la partition entre énergie stockée et énergie dissipée.
Les observations expérimentales montrent que, dans les matériaux métalliques polycristallins, la fraction du travail plastique dissipée en chaleur varie en pratique avec le niveau et la vitesse de déformation (voir par exemple [8]). La plupart de ces mesures expérimentales conduisent à évaluer la partie dissipative de la limite d'écoulement aux alentours de 90% du total. Des deux options exposées ci-dessus, c'est donc la seconde qui semble la plus proche de la réalité physique, en plus d'offrir la meilleure performance numérique.

# 5 Conclusions

Nous avons donc exposé un formalisme variationnel pour le problème incrémental de la plasticité cristalline. Ce formalisme permet d'une part de résoudre les équations constitutives incrémentales sous la forme d'une minimisation sous contraintes, problème pour lequel de nombreux algorithmes d'optimisation sont disponibles. Cette approche permet par ailleurs d'écrire à chaque pas de temps le problème global (problème aux conditions limites) sous forme variationnelle, comme en élasticité.

Nous avons ensuite détaillé différentes options permettant d'intégrer dans ce formalisme des modèles de plasticité cristalline existants. Il est apparu que la solution la plus robuste numériquement est également celle qui est la plus proche des observations expérimentales thermomécaniques. Il serait néanmoins souhaitable de pouvoir disposer d'un modèle constitutif de plasticité cristalline basé directement sur les énergies stockées et dissipées, s'inscrivant dans le cadre variationnel thermomécanique couplé proposé en [11] sans approximations supplémentaires.

# Références

[1] L. Anand, M. Kothari. A computational procedure for rate-independent crystal plasticity. *J. Mech. Phys. Solids*, 44, 525–558, 1996.

[2] E.P. Busso, G. Cailletaud. On the selection of active slip systems in crystal plasticity. *Int. J. Plasticity*, 21, 2212–2231, 2005.

[3] A.M. Cuitiño, M. Ortiz. Computational modeling of single-crystals. *Model. Simul. Mater. Sci. Eng.*, 1, 225–263, 1992.

[4] S.N. Kuchnicki, A.M. Cuitiño, R.A. Radovitzky. Efficient and robust constitutive integrators for single-crystal plasticity modeling. *Int. J. Plasticity*, 22, 1988–2011, 2006.

[5] M. Ortiz, L. Stainier. The variational formulation of viscoplastic constitutive updates. *Comput. Methods Appl. Mech. Engrg*, 171, 419–444, 1999.

[6] D. Peirce, R.J. Asaro, A. Needleman. Material rate dependence and localized deformation in crystalline solids. *Acta Metall.*, 31, 1951–1976, 1983.

[7] J.M. Rice. Inelastic constitutive relations for solids : an internal-variable theory and its application to metal plasticity. *J. Mech. Phys. Solids*, 19, 443–455, 1971.

[8] D. Rittel, A. Bhattacharyya, B. Poon, J. Zhao, G. Ravichandran. Thermomechanical characterization of pure polycrystalline tantalum. *Mater. Sci. Eng. A*, 447, 65–70, 2007.

[9] L. Stainier, A.M. Cuitiño, M. Ortiz. A micromechanical model of hardening, rate sensitivity and thermal softening in bcc single crystals. *J. Mech. Phys. Solids*, 50, 1511–1545, 2002.

[10] G. Weber, L. Anand. Finite deformation constitutive equations and a time integration procedure for isotropic hyperelastic-viscoplastic solids. *Comput. Methods Appl. Mech. Engrg*, 79, 173–202, 1990.

[11] Q. Yang, L. Stainier, and M. Ortiz. A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids. *J. Mech. Phys. Solids*, 54, 401–424, 2006.