# Une famille d'algorithmes robustes pour l'intégration de modèles de plasticité cristalline

Andrei Musienko, Nikolay Osipov, Georges Cailletaud

## ► To cite this version:

Andrei Musienko, Nikolay Osipov, Georges Cailletaud. Une famille d'algorithmes robustes pour l'intégration de modèles de plasticité cristalline. Huitième colloque national en calcul des structures, 2007, Giens, France. pp.271-276. hal-00159693

HAL Id: hal-00159693

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

Submitted on 18 Apr 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 famille d'algorithmes robustes pour
l'intégration de modèles de plasticité
cristalline

Andrey Musienko, Nikolay Osipov, Georges Cailletaud

Centre des Materiaux/Mines Paris, Paristech, CNRS UMR 7633,
B.P. 87, 91003 Evry Cedex, France

RÉSUMÉ. Cet article discute les possibilités d'implémentation numérique des modèles de plasticité cristalline dans un code de calcul par éléments finis. Comme il s'agit d'approches à potentiels multiples, il faut faire face au problème de non-unicité de la solution, puisque plusieurs ensembles de multiplicateurs plastiques (ou pseudo-multiplicateurs en viscoplasticité) sont à même de fournir un incrément de déformation inélastique donné. Plusieurs approches sont testées pour ce qui concerne l'intégration locale et le calcul de la matrice tangente, en petites et en grandes transformations.

ABSTRACT. This paper deals with the numerical implementation of crystal plasticity models into a finite element code. This type of approach involves several potentials, one has then to deal with the problem of non-unicity of the solution, since several sets of plastic multipliers (pseudo-multipliers in viscoplasticity) are able to produce a given inelastic strain rate. Several approaches are tested, both for stress update algorithm, and for the tangent matrix computation, in small and large strain formulation.

MOTS-CLÉS : Algorithme implicite, intégration de loi de comportement, grandes déformations, plasticité cristalline

KEYWORDS: Implicit algorithm, integration of constitutive equations, large transformations, crystal plasticity
# 1. Introduction

The framework of crystal plasticity modeling becomes more and more classical, but there is still a need for more robust and efficient numerical implementation of the material models, due to the specific aspects related to the multiple potentials. The non-unicity of the set of plastic multipliers is well documented in the literature. The problem is related to the fact that crystals often present more than five slip mechanisms, and that several sets of mechanisms can be candidate to accommodate a given inelastic strain rate. A classical solution is provided by viscoplastic models without any threshold (Asaro, 1983, Anand *et al.*, 1994), for which all the slip systems are always active. Time-independent behavior can be recovered by using very high values of the viscosity exponent. Nevertheless, this introduces an artificial regularization near the corners of the yield surface, which can locally influence the direction of the inelastic flow. On the other hand, it must be pointed out that the occurrence of the problem is not related only to the plastic or viscoplastic framework, but to the presence of an elastic region in the stress space. In the present paper, we show a unified view toward this problem, dealing with the problem of unicity for plastic or viscoplastic-with-threshold models, in small or large deformations. In an extended version of the paper, this numerical procedure will be compared to two other classical solutions. The first type of solution comes from the physical analysis of the deformation process, the second one is an algorithmic solution. In the first type of approach, the references are the papers by Bishop and Hill (Bishop *et al.*, 1951), and Chin and Mammel (Chin *et al.*, 1969), who propose a maximization of the plastic work. On the other hand, a series of algorithmic solutions are proposed (Simo *et al.*, 1997, Schröder *et al.*, 1997, Cuitino *et al.*, 1992, Anand *et al.*, 1996, McGinty *et al.*, 2006). A comparative study is provided in (Busso *et al.*, 2005).

A more recent discussion is now open around the dilemma: explicit or implicit integration. Relative performance of implicit versus explicit solvers is tested in (Harewood *et al.*, n.d.). A better computation time is found for implicit solvers, with a factor 1.46 to 5 if compared to explicit ones. Generally speaking, explicit integration is promoted for large scale parallel problems with multiple nonlinear interfaces. An other recent study (Kuchnicki *et al.*, 2006), where the authors have implemented the explicit analogue of the classical implicit method of (Cuitino *et al.*, 1992) concludes that implicit procedures are best suited for quasi-static computations allowing larger time steps. Between the two, one can imagine to perform an explicit Runge-Kutta local integration, together with the use of a consistent tangent matrix (Raphanel *et al.*, 2004). Quadratic convergence is reported by the authors for several specific crystallographic orientations.

One challenge of this work is to review the full implicit solution, in order to clarify (and extend) the practical limits of theoretically expected quadratic convergence. Also, the most simple form of tangent matrix was looked for, in the case of finite strains. Another attempt was made to compare the full implicit and hybrid explicit-implicit solution, to propose the best practical strategy.
## 2. Single crystal model

The model used in this paper was introduced in (Méric et al., 1991), in the small deformation framework. It was extensively used for the computations of polycrystal-line aggregates (Barbe et al., 2001, Diard et al., 2005). Thus only the list of the main equations for the finite strain framework is recalled below. It considers the classical relaxed configuration introduced by Mandel (Mandel, 1973). For the sake of brevity, we do not reproduce here the equations corresponding to kinematic hardening. Classical notations are used.

$$ \mathbf{F} = \underset{\sim}{\mathbf{F}}^{e} \underset{\sim}{\mathbf{F}}^{p} \quad ; \quad \mathbf{E}^{e} = \frac{1}{2}(\underset{\sim}{\mathbf{F}}^{eT} \underset{\sim}{\mathbf{F}}^{e} - \mathbf{I}) \quad ; \quad \mathbf{S} = \underset{\sim}{\mathbf{L}} : \underset{\sim}{\mathbf{E}}^{e} \quad [1] $$

$$ \tau_s = \underset{\sim}{\mathbf{T}} : \underset{\sim}{\mathbf{N}}^s \quad ; \quad \underset{\sim}{\mathbf{T}} = \frac{1}{J} \underset{\sim}{\mathbf{F}}^e \underset{\sim}{\mathbf{S}} \underset{\sim}{\mathbf{F}}^{eT} \quad ; \quad \underset{\sim}{\mathbf{N}}^{s0} = \underset{\sim}{\mathbf{I}}^{s0} \otimes \underset{\sim}{\mathbf{n}}^{s0} \quad [2] $$

$$ \underset{\sim}{\mathbf{L}}^{p} = \underset{\sim}{\dot{\mathbf{F}}}^{p} \underset{\sim}{\mathbf{F}}^{p-1} = \sum_{s} \dot{\gamma}^{s} \underset{\sim}{\mathbf{N}}^{s0} \quad [3] $$

$$ \dot{\gamma}^s = \dot{v}^s \mathrm{sign}(\tau_s) \quad ; \quad \dot{v}^s = \left\langle \frac{|\tau_s| - r_s}{K} \right\rangle^n \quad \text{with } \langle x \rangle = \mathrm{Max}(x, 0) \quad [4] $$

$$ r_s = R_0 + bQ \sum_r h_{sr} q_r = R_0 + Q \sum_r h_{sr} \{1 - e^{-b\nu^r}\} \quad [5] $$

## 3. Integration methods

Two integration methods have been used for this model : an explicit Runge-Kutta method, and an implicit mid-point method resolved by a Newton-Raphson local convergence loop. For both methods, the following set of variables of integration is defined, taking into account the number of slip systems S :

$$ \nu_{\text{int}} = \{ \underset{\sim}{\mathbf{E}}^e, (\nu^s, s = 1..S) \} \quad [6] $$

where the terms $\nu^s$ are the cumulated values of $\gamma_s$. For Runge-Kutta method, one just has to prescribe the rate of each variable. In the implicit case, residuals have to be defined.

For each scheme two types of equations are considered :

– the tensorial equations describing elasto-plastic decomposition. The shape of these equations depends on the chosen framework : for small deformations, (eq. 7) simply results from the elastic-plastic additive partition of the total strain increment; for large deformations, (eq. 8) comes from a first order development of the trial gradient $\underset{\sim}{\mathbf{F}}^* = \underset{\sim}{\mathbf{F}}_{n+1}^e \exp(\Delta\gamma^s \underset{\sim}{\mathbf{N}}^{s0}) \approx \underset{\sim}{\mathbf{F}}_{n+1}^e (\mathbf{I} + \Delta\gamma^s \underset{\sim}{\mathbf{N}}^{s0}) :$

$$ \Delta\underset{\sim}{\mathbf{E}}^e + \sum \Delta\nu^s \mathrm{sign}(\tau_s) \underset{\sim}{\mathbf{n}} = \Delta\underset{\sim}{\mathbf{E}}^{\mathrm{tot}} \quad [7] $$

$$ \Delta\underset{\sim}{\mathbf{E}}^e + \sum \Delta\nu^s \mathrm{sign}(\tau_s) \{\underset{\sim}{\mathbf{N}}^{s0}\underset{\sim}{\mathbf{C}}\} = \Delta\underset{\sim}{\mathbf{E}}^{\mathrm{tot}} \quad [8] $$
– the N scalar equations (N being the number of active slip systems) defining in-elastic flow, which shape changes according to the model (rate-independent (eq.9) or rate-dependent (eq. 10)). For the implicit case one can obtain, for s = 1..N :

$$f^s(\mathbf{E}_{t+\Delta t}^e, (v_{t+\Delta t}^r, r=1..N)) = 0 \quad [9]$$

$$\Delta v^s - \Phi(f^s(\mathbf{E}_{t+\Delta t}^e, (v_{t+\Delta t}^r, r=1..N))) \Delta t = 0 \quad [10]$$

where $f^s$ – criteria of plastic flow and $\Phi(f^s)$ – plastic potential derivative. It has to be checked that all the $\Delta v^s$ are null for systems with a negative value of $f^s$, and strictly positive for systems presenting a viscoplastic ($f^s > 0$) or a plastic ($f^s = 0$) flow.

In the theta-method, the values of all associated forces and parameters are calculated from the internal variables evaluated at an intermediate time $t + \theta\Delta t$, which allows us to introduce the following reduced form of the equations (eq.12), and a Newton-Raphson loop (eq.13) which provides the new estimation of $\Delta v_{int}$ (eq.14) :

$$v_{int}^{\theta} = v_{int}^{t} + \theta \Delta v_{int} \quad 0 \le \theta \le 1 \quad [11]$$

$$\mathcal{F}(v_{int}^{\theta}, \Delta v_{int}) = \mathcal{F}_0 \quad [12]$$

$$\mathcal{R}^k = \mathcal{F}^k - \mathcal{F}_0 \quad [13]$$

$$\Delta v_{int}^{k+1} = \Delta v_{int}^{k} - [\Delta \mathcal{F}_{\theta}^{k+1}]^{-1} \mathcal{R}^{k} \quad [14]$$

where the $\theta$ subscript indicates calculation of $\Delta v_{int}$ at time $t + \theta\Delta t$. In fact, the method provides the most robust integration schemes for $\theta = 1$. After convergence, the top left subcomponent of the inverted Jacobian matrix $[\Delta \mathcal{F}_{\theta}^{k+1}]^{-1}$ relates the change in $\Delta \mathbf{E}_{\sim}^e$ with respect to a change in $\Delta \mathbf{E}_{\sim}^{tot}$. The elastic constitutive equations can then be used to evaluate the partial derivatives of $\Delta \sigma$ with respect to $\Delta \mathbf{E}_{\sim}^{tot}$. Therefore, the non-linear tangent matrix consistent with the integration algorithm is (Simo et al., 1997) :

$$\mathbf{D}_{\sim} = \frac{\delta\Delta\sigma_{\sim}}{\delta\Delta\mathbf{E}_{\sim}^{tot}} = \mathbf{C}_{\sim} [\Delta\mathcal{F}_{\theta}^{k+1}]_{e}^{-1} \quad [15]$$

An alternative to the full implicit algorithm can be considered by using and *hybrid* integration technique. In such a case, the local integration is performed by a Runge-Kutta method, but a tangent matrix is then computed. The construction of this matrix is allowed by a "one-shot" resolution of the implicit system built with the residuals attached to the theta-method, which is initialized by means of the solution of the Runge-Kutta integration. In this case, the tangent matrix is not exactly consistent with the integration scheme, but its quality is still quite satisfactory and usually provides quadratic convergence.
Let us mention that a simplified form of tangent matrix has also be successfully
tested in our large deformation computations :

$$
\boldsymbol{\sigma}^{\nabla} = \left( \frac{1}{J} \boldsymbol{\phi}_{fw}^{\sim e} \left[ \boldsymbol{L}^{\sim} : \left( \frac{\partial \boldsymbol{f}_e}{\partial \Delta \boldsymbol{E}_{\sim e}} \right)^{-1} \right] - \boldsymbol{M}_{L}^{\sim} (\boldsymbol{\sigma}_{\sim}) - \boldsymbol{M}_{R}^{\sim} (\boldsymbol{\sigma}_{\sim}) \right) : \boldsymbol{D}_{\sim} \quad [16]
$$

4. Discussion and concluding remarks

It is worth noting that the proposed approach does not need any specific technique
for the choice of the slip systems. This produces faster convergences, if compared to
solutions studied by other authors ((Anand et al., 1996) with a SVD decomposition of
the matrix; (McGinty et al., 2006) with an iterative scheme that introduces the active
slip systems one-by-one). The method is more inspired from (Simo et al., 1997) who
proposes just a control of the active slip systems at each iteration. All the techniques
were compared together in a previous paper, and were found to give similar results,
even for non proportional loading paths (Busso et al., 2005). It has to be noted that the
literature shows very few numerical tests with arbitrary crystal orientations : authors
use to try either single slip case, or perfect multiple slip. In order to investigate the way
the method performs in a polycrystal, a series of randomly chosen crystal orientations
have been chosen, in a unique cube element. A simple interaction matrix [h] is chosen,
with h<sub>ij</sub> = h<sub>2</sub> + (1 - h<sub>2</sub>)δ<sub>ij</sub>. Different values are tested for h<sub>2</sub> (resp. 0, 0.5, 1, 1.5). A
tension from 0 to 1% is applied to the cube. The global convergence is illustrated in
Fig.4. It is found that the convergence is quadratic, as shown in Fig.4a for most of
the cases (small dots in the standard triangle, Fig.4c). Nevertheless, the convergence
may become linear only for orientations near multiple slip (large dots in the standard
triangle, and Fig.4b), since, in this case, the number of active slip systems may vary
during the global equilibrium iterations.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Examples of quadratic (a) and linear (b) convergence according to (c) the position in the standard triangle – inverse pole figure (tension direction) with 40 random orientations. Bold points – linear global convergence, small dots – quadratic convergence</figcaption></figure>
**5. Bibliographie**

Anand L., Kalidindi S., « The process of shear band formation in plane strain compression of FCC metals : Effects of crystallographic texture », *Mech. of Materials*, vol. 17, p. 223-243, 1994.

Anand L., Kothari M., « A computational procedure for rate-independent crystal plasticity », *J. Mech. Phys. Sol.*, vol. 44, n° 4, p. 525-558, 1996.

Asaro R., « Crystal Plasticity », *J. of Applied Mechanics*, vol. 50, p. 921-934, 1983.

Barbe F., Forest S., Cailletaud G., « Intergranular and intragranular behavior of polycrystalline aggregates. Part 2 : Results », *Int. J. of Plasticity*, vol. 17, n° 4, p. 537-563, 2001.

Bishop J., Hill R., « A Theoretical Derivation of the Plastic Properties of a Polycrystalline Face-Centered Metal », *Philosophical Magazine*, vol. 42, p. 414-427, 1951.

Busso E., Cailletaud G., « On the selection of active slip systems in crystal plasticity », *Int. J. of Plasticity*, 2005.

Chin G., Mammel W., « Generalization and Equivalence of the Minimum Work (Taylor) and Maximum Work (Bishop-Hill) Principles for Crystal Plasticity », *Trans. of the Met. Soc. of AIME*, vol. 245, p. 1211-1214, 1969.

Cuitino A., Ortiz M., « Computational modelling of single crystals », *Modelling Simul. Mater. Sci. Eng.*, vol. 1, p. 225-263, 1992.

Diard O., Leclercq S., Rousselier G., Cailletaud G., « Evaluation of finite element based analysis of 3D multicrystalline aggregates plasticity : Application to crystal plasticity model identification and the study of stress and strain fields near grain boundaries », *Int. J. of Plasticity*, vol. 21, n° 4, p. 691-722, 2005.

Harewood F. J., McHugh P. E., « Comparison of the implicit and explicit finite element methods using crystal plasticity », *Computational Materials Science*, n.d.

Kuchnicki S. N., Cuitino A. M., Radovitzky R. A., « Efficient and robust constitutive integrators for single-crystal plasticity modeling », *International Journal of Plasticity*, vol. 22, n° 10, p. 1988-2011, 2006.

Mandel J., « Equations constitutives et directeurs dans les milieux plastiques et viscoplastiques », *Int. J. Solids Structures*, vol. 9, n° 6, p. 725-740, 1973.

McGinty R. D., McDowell D. L., « Semi-implicit integration scheme for rate independent finite crystal plasticity », *International Journal of Plasticity*, vol. 22, n° 6, p. 996-1025, 2006.

Méric L., Cailletaud G., « Single Crystal Modeling for Structural Calculations. Part 2 : Finite Element Implementation », *J. of Engng. Mat. Technol.*, vol. 113, p. 171-182, 1991.

Raphanel J. L., Ravichandran G., Leroy Y. M., « Three-dimensional rate-dependent crystal plasticity based on Runge-Kutta algorithms for update and consistent linearization », *International Journal of Solids and Structures*, vol. 41, n° 22-23, p. 5995-6021, 2004.

Schröder J., Miehe C., « Aspects of computational rate-independent crystal plasticity », *Computational Materials Science*, vol. 9, p. 168-176, 1997.

Simo J., Hughes T., *Computational Inelasticity*, Springer Verlag, 1997.