# A generic formulation of anisotropic thermo-elastoviscoplasticity at finite deformations for finite element codes

Mohamed Abatour<fnref n="1" />,<fnref n="2" /> ID · Kais Ammar<fnref n="1" /> · Samuel Forest<fnref n="1" /> ID · Cristian Ovalle<fnref n="1" /> ID · Nikolay Osipov<fnref n="2" /> · Stéphane Quilici<fnref n="2" />

Received: 21 January 2024 / Accepted: 15 August 2024
© The Author(s) 2024

## Abstract

Most commercial Finite Element codes rely on hypo-elastoviscoplastic constitutive equations for structural computations at large deformations. Such formulations are known to suffer from physical and thermodynamical flaws but their versatility and ease of implementation made them ubiquitous. The paper presents an alternative systematic and thermodynamically consistent extension of anisotropic thermo-elastoviscoplastic constitutive equations at finite strain. The formulation is based on the well-known multiplicative decomposition of the deformation gradient into a thermoelastic and an inelastic part. This decomposition relies on Mandel's isoclinic intermediate configuration. The present framework covers a wide range of multi-mechanism elastoplastic models so that it can advantageously replace the constitutive part of existing codes. The choice of a suitable hyperelastic potential, hardening variables and anisotropy evolution laws are discussed. The concept of plastic spin is used and can be either derived from general representation theorems or obtained from the normality rule. The effect of the plastic spin is discussed in the light of examples involving anisotropic plasticity. The response of the proposed formulation is compared to hypoelastic models in the case of several structural applications. The implementation of this methodology in a commercial FE object-oriented code is detailed. We show how to extend readily a wide range of small strain nonlinear constitutive models to finite deformations. The paper contains original features such as multimechanism based inelastic contributions, a Mandel-stress tensor based nonlinear kinematic hardening rule containing dynamic and static recovery terms, and a study of lattice rotation in single crystalline volume elements and turbine blades. The implementation is proved to be competitive with respect to existing hypo-elastoviscoplastic formulations in terms of CPU time.

**Keywords** Finite deformations · Finite element · Elastoplasticity · Plastic spin · Anisotropic plasticity

This article is dedicated to Jörg Schröder on the occasion of his 60th birthday. The subject of the article reflects 30 years of Franco-German complicity in the field of computational mechanics and constitutive laws for materials undergoing finite deformations.

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

Samuel Forest
samuel.forest@minesparis.psl.eu;
samuel.forest@mines-paristech.fr

Mohamed Abatour
mohamed.abatour@transvalor.com

Kais Ammar
kais.ammar@minesparis.psl.eu

Cristian Ovalle
cristian.ovalle@minesparis.psl.eu

Nikolay Osipov
nikolay.osipov@transvalor.com

Stéphane Quilici
stephane.quilici@transvalor.com

## 1 Introduction

Finite strain elastoplasticity, as observed in various materials, requires combined geometric and material nonlinear analysis of solids. Since the early sixties, a myriad of strategies have been developed to extend the well-established infinitesimal elastoplasticity theory to finite transformations [1, 2]. A widely used approach is the so-called hypoelastic formulation relying on additive decomposition of the total deformation rate into elastic and inelastic parts, and constitutive equations for objective stress rates [3, 4]. However, this framework has been the subject of much controversy.

<footnote n="1">¹ Mines Paris, Centre des matériaux, CNRS UMR 7633, PSL University, 91003 Évry, France</footnote>

<footnote n="2">² TRANSVALOR S.A., Centre d'affaires la Boursidière, 92357 Le Plessis-Robinson, France</footnote>
First, the constitutive equations are generally not integrable, which results in spurious energy dissipation in the elastic regime even prior to yield (see, e.g. [5–7] among others). In order to recover the integrability of hypoelastic formulations, a logarithmic rate has been put forward by [8, 9]. However, it was shown in [10] that the post-yield response of the logarithmic rate is inconsistent with the notion of elasticity. Second, the rate of elastic deformation is related to a non unique objective stress rate [11, 12]. In addition, some hypoelastic models are well-known to depend on the reference configuration, e.g. Green–Naghdi, and logarithmic rates [13]. Formulations relying upon the additive decomposition of the Green-Lagrange strain tensor were developed [1]. Given that the considered strain measure is symmetric, the latter model cannot describe full anisotropy. Another formulation based on the additive split of logarithmic strain was suggested by [14]. These formulations however imply a dependence on the choice of the reference configuration, meaning that the form of the constitutive laws is not left unchanged by the change of reference configuration [13]. Significant differences in the plastic strain localization behaviour of material laws involving additive and multiplicative decompositions were recently evidenced in [15].

Since the aforementioned models are unsatisfactory in describing finite deformations, several authors have put forward the necessity to introduce the multiplicative decomposition as a general way of describing the kinematics at finite strain [16, 17]. This decomposition assumes the existence of an intermediate configuration that is generally not unique. The concept of isoclinic intermediate configuration was proposed first by [18, 19]. Since then, several models relying upon the multiplicative decomposition and a hyperelastic potential for the stress have been developed (e.g. [7, 20–23]). It is noteworthy to mention that there exist different versions of the multiplicative decomposition, e.g. considering a decomposition in the reverse order [24, 25], decomposition assuming symmetric elastic part or symmetric plastic part [17, 26]. Alternatively, a rigorous presentation of this theory has been proposed by [27] based on material isomorphisms. This approach defines a special class of materials described by elastic ranges that are independent of plastic deformation.

Constitutive modeling of kinematic hardening is still an active area of research even in the small strain regime [28]. Kinematic hardening models at finite strain can be generally classified into two main groups; both are regarded as extensions of the Armstrong-Frederick model [29]. The first one employs an evolution equation for the back stress, or the so-called 'Chaboche-type' model [30]. The second one involves an additional multiplicative decomposition of the plastic part of the deformation gradient into storage and dissipative parts [20, 23, 31–33]. The difference between these models becomes visible, particularly in cases where principal axes rotate, e.g. simple glide and torsion [34]. Significant dif-

ferences are observed for the special case of linear kinematic hardening (Prager model) which leads to stress oscillations for Jaumann rate [34]. Meanwhile, many studies have shown that different extensions of the Armstrong-Frederick model yield, at least qualitatively, similar results [32, 35, 36]. The differences are mainly due to second-order effects for loading conditions involving large rotations, apparent in simple shear [32, 37].

At finite strain, the material does not undergo only macroscopic stretches and rotations but also rotations of its substructure described by some privileged directions called *directors*. The notion of plastic spin describes the evolution of material's directors with plasticity. As highlighted by several works (e.g. [38, 39]), the plastic spin is undetermined for isotropic materials and is often assumed to vanish. Several constitutive equations for the plastic spin have been proposed independently by [40, 41]. It has been argued that multiple plastic spins are required since each internal variable has a different nature and, consequently, a different rotating frame is required for each internal variable [42]. Furthermore, the plastic spin is useful to prevent stress oscillation, particularly in the case of the simple shear problem [41] when applied to elastoplastic materials with linear kinematic hardening. Constitutive modeling of plastic spin can be achieved through: (i) an additional ad-hoc constitutive equation [40, 43, 44] (ii) generalized normality conditions [19, 45, 46].

The additive hypo-elastoplastic models are widely used in finite element codes for the sake of computational ease. Small deformation material models are adapted to large deformation problems with appropriate tangent operators and stress/strain measures. To our knowledge, there is no finite element software supporting a systematic extension of infinitesimal models using the multiplicative decomposition. Meanwhile, this approach is already used, exclusively, for some particular models, e.g. in MSC Marc [47], or by using user material subroutines in ABAQUS [7, 23].

This work aims to propose a generic and systematic extension of small strain models to finite deformations. The multiplicative decomposition of the deformation gradient into thermoelastic and inelastic parts is adopted. All constitutive equations are expressed in the local isoclinic intermediate configuration. Accordingly, the stress is related to the elastic strain through a hyperelastic law. Emphasis is put on the broad range of addressed constitutive equations and on the numerical implementation of the present formulation in a commercial finite element code provided with an object-oriented interface [48]. By this approach, small strain models, including isotropic and anisotropic plasticity models with various isotropic and kinematic hardening rules, can be reused in the framework of the multiplicative decomposition. On the other hand, unlike usual hypoelastic models, various constitutive equations for the plastic spin can be formulated within this framework. A significant part of the paper is a
review of existing formulations and a systematic comparison of the results obtained with the proposed formulation to previous ones. Original features deal with the formulation of kinematic hardening with dynamic and static recovery terms, coupling of various inelastic mechanisms, and detailed description of the rotation of directors in single and polycrystals.

The manuscript is organized as follows. Section 2 presents a general thermodynamical framework for constitutive modeling of anisotropic finite strain thermo-elastoviscoplasticity based on the multiplicative decomposition. Kinematic assumptions and a thermodynamically consistent derivation of constitutive equations for multi-mechanism modeling are depicted. Several models accounting for kinematic hardening are discussed. Two constitutive choices of plastic spin are presented. In Sect. 3, we describe the implementation of the present formulation in a commercial finite element code by taking advantage of its object oriented interface. Special attention is paid to the integration of constitutive equations and the construction of tangent matrices. Two integration schemes are adopted: using the exponential mapping or correcting residuals in order to fulfill the plastic incompressibility condition. Finally, the response of some particular models is presented in Sects. 4 and 5 for a volume element and structural applications, respectively. Through these applications, comparisons are made with hypoelastic formulations presented in the Appendix A. For the sake of brevity, the thermomechanical coupling of the balance laws is not considered in the present work. This means that the temperature fields are assumed to be given in all samples, for instance after solving independently the heat equation, as often done in engineering computations.

The notations used throughout are as follows. $\mathbf{a}$, $\mathbf{A}$, and $\mathbf{A}^T$ stand respectively for first, second and fourth order tensors. The transpose, inverse, transpose of inverse and material time derivative of a second order tensor are denoted by $\tilde{\mathbf{A}}^T$, $\tilde{\mathbf{A}}^{-1}$, $\tilde{\mathbf{A}}^{-T}$, and $\tilde{\mathbf{A}}$ respectively. Double contractions are denoted by $\tilde{\mathbf{A}} : \tilde{\mathbf{B}} = A_{ij}B_{ij}$ and $\tilde{\mathbf{A}} : \tilde{\mathbf{B}} = A_{ijkl}B_{kl} \mathbf{e}_i \otimes \mathbf{e}_j$, with $(\mathbf{e}_i)_{i=1,2,3}$ being a Cartesian orthonormal basis. The following tensor products are used: $\tilde{\mathbf{A}} \otimes \tilde{\mathbf{B}} = A_{ij}B_{kl} \mathbf{e}_i \otimes \mathbf{e}_j \otimes \mathbf{e}_k \otimes \mathbf{e}_l$, $\tilde{\mathbf{A}} \otimes \tilde{\mathbf{B}} = A_{ik}B_{jl} \mathbf{e}_i \otimes \mathbf{e}_j \otimes \mathbf{e}_k \otimes \mathbf{e}_l$, $\tilde{\mathbf{A}} \otimes \tilde{\mathbf{B}} = A_{il}B_{jk} \mathbf{e}_i \otimes \mathbf{e}_j \otimes \mathbf{e}_k \otimes \mathbf{e}_l$.

# 2 General framework for elastoviscoplastic modeling at finite strains

## 2.1 Kinematics

The deformation gradient is multiplicatively split as

$$ \tilde{\mathbf{F}} = \tilde{\mathbf{F}}^e \tilde{\mathbf{F}}^p \qquad (1) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Illustration of the local isoclinic intermediate configuration</figcaption></figure>

where $\tilde{\mathbf{F}}^e$ and $\tilde{\mathbf{F}}^p$ denote the thermoelastic and the inelastic parts, respectively, see Fig. 1. In contrast to $\tilde{\mathbf{F}}$, the parts $\tilde{\mathbf{F}}^e$ and $\tilde{\mathbf{F}}^p$ are not necessarily defined as gradients of one-to-one mappings. As a consequence of possible incompatibility of thermoelastic and plastic deformation fields, the intermediate configuration is local, i.e. it belongs to the immediate vicinity of a material point and is obtained from the current deformed configuration by a purely elastic unloading of this vicinity only. The volume changes due to elastic and plastic deformations are respectively denoted by

$$ J_e = \det(\tilde{\mathbf{F}}^e) = \frac{\rho_e}{\rho}, \quad J_p = \det(\tilde{\mathbf{F}}^p) = \frac{\rho_0}{\rho_e} \qquad (2) $$

and $\rho$, $\rho_e$ and $\rho_0$ stand for the mass densities in current, intermediate and reference configurations, respectively. The decomposition (1), however, is not unique since any invertible transformation $\tilde{\mathbf{H}}$ can be introduced such that

$$ \tilde{\mathbf{F}} = (\tilde{\mathbf{F}}^e \tilde{\mathbf{H}})(\tilde{\mathbf{H}}^{-1} \tilde{\mathbf{F}}^p) = \tilde{\mathbf{F}}^{e*} \tilde{\mathbf{F}}^{p*} \qquad (3) $$

To remedy this concern, we shall define some material directors. These privileged directors are attached to some microstructural features and undergo a different transformation than the material. For solids, we can always define a triad of directors describing its microstructure (fibers, crystal lattice vectors,...). In general, an infinity of directors may be identified for a material. Nevertheless, if the relative spin of all directors with respect to a single triad of directors is known, the problem shall be simplified by considering only this triad. It leads to the concept of local isoclinic intermediate configuration, advocated in [19], where the material directors have the same inclination or orientation as in the reference configuration. The intermediate configuration is obtained by elastic virtual unloading, whereas the rigid-body rotation of directors is included in the thermoelastic part of the deformation gradient. The isoclinic local configuration is uniquely
defined up to a symmetry operation belonging to the symmetry group of the material. It is one example of the co-rotated intermediate configurations recommended in [49] where the used structural tensors are also based on a triad of directors [50].

The velocity gradient is additively decomposed in the current configuration

$$ \dot{\underline{F}} \underline{\underline{F}}^{-1} = \underline{\underline{L}}^e + \underline{\underline{F}}^e \underline{\underline{L}}^p \underline{\underline{F}}^{e-1} \quad (4) $$

$\underline{\underline{L}}^e$ is the purely elastic part of the velocity gradient and $\underline{\underline{L}}^p$ is referred to as the plastic deformation rate in the intermediate configuration. The latter, in turn, can be split into symmetric and skew-symmetric parts as

$$ \begin{cases} \underline{\underline{L}}^e = \dot{\underline{\underline{F}}}^e \underline{\underline{\underline{F}}}^{e-1} = \underline{\underline{D}}^e + \underline{\underline{W}}^e \\ \underline{\underline{L}}^p = \dot{\underline{\underline{F}}}^p \underline{\underline{\underline{F}}}^{p-1} = \underline{\underline{D}}^p + \underline{\underline{W}}^p \end{cases} \quad (5) $$

where $\underline{\underline{D}}^e = \text{sym}(\underline{\underline{L}}^e)$, $\underline{\underline{D}}^p = \text{sym}(\underline{\underline{L}}^p)$ are the elastic and the plastic strain rates, and $\underline{\underline{W}}^e = \text{skw}(\underline{\underline{L}}^e)$ and $\underline{\underline{W}}^p = \text{skw}(\underline{\underline{L}}^p)$ are the elastic and plastic spin tensors.

## 2.2 Thermodynamic framework

The local form of the Clausius-Duhem inequality expressed in the local current configuration is written as

$$ \frac{\underline{\underline{\sigma}} : \underline{\underline{D}}}{\rho} - (\dot{\underline{\psi}} + \dot{T}\eta) - \frac{1}{\rho} \frac{\underline{\underline{q}} \cdot \underline{\underline{g}}}{T} \geq 0, \quad \underline{\underline{g}} = \nabla_x T \quad (6) $$

where $\underline{\underline{q}}$ is the heat flux vector. This inequality involves the Helmholtz free energy and entropy densities per unit mass. The former is defined as $\underline{\underline{\psi}} = e - T\eta$ where $e$ is the specific internal energy function, $T$ is the absolute temperature and $g$ is the current spatial gradient of temperature. The volume density of internal forces w.r.t. the intermediate configuration is given by

$$ J_e \underline{\underline{\sigma}} : \underline{\underline{D}} = \underline{\underline{\Pi}}^e : \dot{\underline{\underline{E}}}^e + \underline{\underline{M}} : \underline{\underline{L}}^p \quad (7) $$

where

$$ \begin{cases} \underline{\underline{E}}^e = \frac{1}{2} (\underline{\underline{F}}^{eT} \underline{\underline{F}}^e - \underline{\underline{1}}) \\ \dot{\underline{\underline{E}}}^e = \underline{\underline{\underline{F}}}^{eT} \underline{\underline{\underline{D}}}^e \underline{\underline{\underline{F}}}^e \\ \underline{\underline{\underline{\Pi}}}^e = J_e \underline{\underline{\underline{F}}}^{e-1} \underline{\underline{\underline{\sigma}}} \underline{\underline{\underline{F}}}^{e-T} \\ \underline{\underline{\underline{M}}} = J_e \underline{\underline{\underline{F}}}^{eT} \underline{\underline{\underline{\sigma}}} \underline{\underline{\underline{F}}}^{e-T} = \underline{\underline{C}}^e \underline{\underline{\underline{\Pi}}}^e \end{cases} \quad (8) $$

$\underline{\underline{\Pi}}^e$ and $\underline{\underline{M}}$ denote the Piola and Mandel stress tensors, expressed in the intermediate configuration. The Mandel stress tensor $\underline{\underline{M}}$ is generally non-symmetric, in contrast to

$\underline{\underline{\Pi}}^e$. By expressing the dissipation inequality in the intermediate configuration, we obtain

$$ \underline{\underline{\Pi}}^e : \dot{\underline{\underline{E}}}^e + \underline{\underline{M}} : \underline{\underline{L}}^p - \rho_e (\dot{\underline{\psi}} + T\eta) - \frac{\underline{\underline{q}}_e \cdot \underline{\underline{g}}_e}{T} \geq 0 \quad (9) $$

where $\underline{\underline{q}}_e = J_e \underline{\underline{\underline{F}}}^{e-1} \underline{\underline{q}}$ and $\underline{\underline{g}}_e = \underline{\underline{\underline{F}}}^{eT} \underline{\underline{g}}$. The specific free energy density $\underline{\underline{\psi}}(\underline{\underline{E}}^e, T, \underline{\underline{\alpha}}_I)$ is a function of elastic strain $\underline{\underline{E}}^e = (\underline{\underline{\underline{F}}}^{eT} \underline{\underline{\underline{F}}}^e - \underline{\underline{1}})/2$, i.e. the Green-Lagrange strain measure with respect to the intermediate configuration, temperature $T$ and internal variables $\underline{\underline{\alpha}}_I$ which are scalar and/or tensor quantities accounting for hardening properties. It follows that

$$ \dot{\underline{\psi}} = \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{E}}^e} : \dot{\underline{\underline{E}}}^e + \sum_I \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{\alpha}}_I} \dot{\underline{\alpha}}_I + \frac{\partial \underline{\underline{\psi}}}{\partial T} \dot{T} \quad (10) $$

The Clausius-Duhem inequality takes the form:

$$ \left( \underline{\underline{\Pi}}^e - \rho_e \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{E}}^e} \right) : \dot{\underline{\underline{E}}}^e + \underline{\underline{M}} : \underline{\underline{L}}^p - \rho_e \sum_I \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{\alpha}}_I} \dot{\underline{\alpha}}_I - \rho_e \left( \frac{\partial \underline{\underline{\psi}}}{\partial T} + \eta \right) \dot{T} - \frac{\underline{\underline{q}}_e \cdot \underline{\underline{g}}_e}{T} \geq 0 \quad (11) $$

The following state laws are adopted

$$ \underline{\underline{\Pi}}^e = \rho_e \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{E}}^e}, \quad \eta = -\frac{\partial \underline{\underline{\psi}}}{\partial T} \quad (12) $$

so that the intrinsic dissipation remains as

$$ \phi_{in} = \underline{\underline{M}} : \underline{\underline{L}}^p - \rho_e \sum_I \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{\alpha}}_I} \dot{\underline{\alpha}}_I \geq 0 \quad (13) $$

together with the condition $\underline{\underline{q}}_e \cdot \underline{\underline{g}}_e \leq 0$. The condition (13) is satisfied for any process if there exists a convex potential $\Omega$ w.r.t. the generally non-symmetric second order tensor $\underline{\underline{M}}$ and concave w.r.t. $\underline{\underline{A}}_I$ such that

$$ \underline{\underline{L}}^p = \frac{\partial \Omega(\underline{\underline{M}}, \underline{\underline{A}}_I)}{\partial \underline{\underline{M}}} , \quad \dot{\underline{\alpha}}_I = - \frac{\partial \Omega(\underline{\underline{M}}, \underline{\underline{A}}_I)}{\partial \underline{\underline{A}}_I} \quad (14) $$

$\underline{\underline{A}}_I$ are the thermodynamic forces associated with the state variables $\underline{\underline{\alpha}}_I$:

$$ \underline{\underline{A}}_I = \rho_e \frac{\partial \underline{\underline{\psi}}}{\partial \underline{\underline{\alpha}}_I} \quad (15) $$

In the case of time-independent plasticity, the flow rule in Eq. (14) is rewritten

$$ \underline{\underline{L}}^p = \dot{\lambda} \frac{\partial f(\underline{\underline{M}}, \underline{\underline{A}}_I)}{\partial \underline{\underline{M}}} \quad (16) $$
where $f(M, A_I)$ is the yield function and $\dot{\lambda}$ denotes the plastic multiplier which can be determined by use of the consistency condition as

$$ \dot{f}(M, A_I) = \frac{\partial f}{\partial \tilde{M}} : \dot{\tilde{M}} + \frac{\partial f}{\partial A_I} : \dot{A}_I = 0 \quad (17) $$

For viscoplasticity with a threshold, from Eqs. (16) and Eq. (14), a viscoplastic multiplier can be defined as

$$ \dot{\lambda} = \frac{\partial \Omega(M, A_I)}{\partial f(M, A_I)} \quad (18) $$

The existence of a convex potential from which the flow rule and the evolution law of internal variables are derived is sufficient to satisfy the dissipation inequality in Eq. (13), but it is not necessary.

## 2.3 Thermo-hyperelasticity

A myriad of free energy potentials are available to model hyperelasticity at finite strains. A widely used free energy potential is the so-called St. Venant–Kirchhoff model, as a straightforward generalization of Hooke’s law

$$ \rho_e \psi^e(\tilde{E}^e) = \frac{1}{2} \tilde{E}^e : \tilde{C} : \tilde{E}^e, \quad \text{with} \quad \tilde{E}^e = \frac{1}{2}(\tilde{F}^e T \tilde{F}^e - \tilde{1}) \quad (19) $$

where $\tilde{C}$ is the fourth-order elasticity moduli. This non-polyconvex potential [51] fails to respond appropriately in some cases e.g. the stress needed to shrink a bar to zero volume goes to zero which is physically unreasonable [52]. Several modified versions of this model are proposed in literature in order to circumvent the aforementioned limitation [52–54] e.g. the isotropic neo-Hookean model

$$ \rho_e \psi^e(\tilde{E}^e) = \frac{1}{2} \lambda(\log(J_e))^2 + \mu (\text{trace}(\tilde{E}^e) - \log(J_e)) \quad (20) $$

$\lambda$ and $\mu$ are Lamé coefficients. We mention that the proposed formulation is not restricted to a unique hyperelastic model. The model in Eq. (19) is proposed by default. In particular this choice is sufficient to represent the elastic part of the elastoplastic behavior of metals and alloys which are characterized by small elastic strains.

Two approaches for introducing thermoelasticity at finite strain may be distinguished. In the first approach, two configurations of material sample are considered: The initial configuration at uniform reference temperature, and the deformed configuration characterized by non-uniform stress and temperature fields [27, 55, 56]. The second approach has been proposed by [55, 57–59] and by [20, 60, 61] in the

framework of thermo-elastoplasticity, where a supplementary intermediate configuration is considered. In the purely thermoelastic case, the deformation gradient is then split multiplicatively into thermal and elastic parts as

$$ \tilde{F} = \tilde{F}^{el} \tilde{F}^{\theta} \quad (21) $$

An alternative decomposition in the form of $\tilde{F} = \tilde{F}^{\theta} \tilde{F}^{el}$ is suggested by [62]. For isotropic materials, the two approaches yield identical or similar results [59]. When extended to anisotropic elastoplasticity, the latter approach gives rise to potential decompositions of the form $\tilde{F} = \tilde{F}^{el} \tilde{F}^{\theta} \tilde{F}^{p}$ and all permutations of this decomposition can be found in the literature<fnref n="1" /> without unambiguous justification for the best-suited choice. That is why the first approach, see [56], which relies on a unified thermoelastic deformation $\tilde{F}^e$ is preferred in the following, so that superfluous sub-decompositions are avoided.

The thermoelastic part of the Helmholtz free energy $\psi^{the}$ is defined as a function of the thermoelastic deformation tensor, still defined by Eq. (19), and temperature as

$$ \rho_e \psi^{the}(\tilde{E}^e, T) = \rho_e \psi^e(\tilde{E}^e) - \Delta T \beta : \tilde{E}^e - \frac{\rho_e C_\varepsilon}{2T_0} \Delta T^2 \quad (22) $$

where $\beta$ is a second order symmetric stress temperature tensor [63] and $C_\varepsilon$ is the specific heat at constant strain. The difference $\Delta T = T - T_0$ depends on a constant reference parameter $T_0$. If the free energy definition (19) is adopted, the stress is obtained as

$$ \tilde{\Pi}^e = \tilde{C} : (\tilde{E}^e - \tilde{\alpha}_T \Delta T), \quad \tilde{\alpha}_T = \tilde{C}^{-1} : \tilde{\beta} \quad (23) $$

where $\tilde{\alpha}_T$ denotes the thermal expansion second order tensor which, for isotropic materials, can be reduced to only one scalar parameter as $\tilde{\alpha}_T = \alpha_T \tilde{1}$. Due to the fact that in usual materials thermal expansion always remains small, the proposed framework closely related to the usual small strain formulation will be sufficient.

## 2.4 Multimechanism dissipation potential

The inelastic deformation observed on the macroscale has various origins at the level of material microstructure (dislocation slip, twinning, grain boundary sliding, viscosity and molecular orientation in polymers...). Each mechanism leads to a specific type of non-linearity and constitutive equations. To sum up the contributions of all mechanisms, the total dis-

<footnote n="1">¹ $\tilde{F} = \tilde{F}^{el} \tilde{F}^{p} \tilde{F}^{\theta}$ according to [20] and $\tilde{F} = \tilde{F}^{el} \tilde{F}^{\theta} \tilde{F}^{p}$ in [60, 61]</footnote>
sipation potential is written in the form

$$ \Omega(\tilde{\mathbf{M}}, \mathbf{A}_I) = \sum_{k}^{N} \Omega_k(\tilde{\mathbf{M}}, \mathbf{A}_I) \quad (24) $$

$\Omega_k$ denotes the potential associated with the individual viscoplastic mechanism $k$. Each mechanism is characterized by several internal variables summarizing, at a given time, the material state and the influence of the past thermomechanical loading. The yield surface is defined in the stress-hardening variable (and temperature) space. For a given temperature and hardening, the elastic range is limited by the yield surface. It is a part of the vector space of dimension 9 of non-symmetric second order tensors denoted by $\mathcal{D}_e = \{\tilde{\mathbf{M}}/f(\tilde{\mathbf{M}}, \mathbf{A}_I, T) \le 0\}$. The condition $f(\tilde{\mathbf{M}}, \mathbf{A}_I, T) = 0$ defines the yield surface and is chosen here of the form

$$ f(\tilde{\mathbf{M}}, \tilde{\mathbf{X}}, R) = [\tilde{\mathbf{M}} - \tilde{\mathbf{X}}]_{eq} - R_0 - R(p, \alpha_I) \quad (25) $$

where $[\diamond]_{eq}$ denotes an equivalent stress measure involving appropriate invariants of the tensor inside the brackets, $R_0$ is the initial yield stress, $R$ describes the isotropic hardening law depending on the accumulated plastic strain $p$ and internal variables $\alpha_I$, and $\tilde{\mathbf{X}}$ is the back stress. Considering for instance the case of two internal variables $\alpha_I = (\alpha, r)$ and using Eq. (15), the associated forces $\tilde{\mathbf{X}}$ and $R$ are obtained through

$$ \tilde{\mathbf{X}} = \rho_e \frac{\partial \psi}{\partial \tilde{\alpha}}, \quad R = \rho_e \frac{\partial \psi}{\partial r} \quad (26) $$

where $\tilde{\alpha}$ and $r$ are the internal variables associated to kinematic and isotropic hardening, respectively. Multiple yield functions and multiple kinematic hardening variables can be used in the context of multimechanism approach, as proposed in [64, 65]

It follows from Eqs. (16) and (18) that the inelastic strain rate is the sum of the individual contributions of all mechanisms:

$$ \bar{\mathbf{L}}^p = \sum_k \frac{\partial \Omega_k(\tilde{\mathbf{M}}, \mathbf{A}_I)}{\partial \tilde{\mathbf{M}}} \quad (27) $$

This method circumvents the decomposition of the inelastic deformation into various contributions, often used in the literature in the form $\tilde{\mathbf{F}}^p = \tilde{\mathbf{F}}^p l \tilde{\mathbf{F}}^v \tilde{\mathbf{F}}^d \tilde{\mathbf{F}}^{tr}...$ including plastic (rate-independent), viscoplastic, damage or transformation deformations. Such decompositions and all their possible permutations are hard to justify or define unambiguously. In contrast, each mechanism contributes incrementally to the inelastic deformation rate and can generally not be time-integrated into one single deformation part.

## 2.5 Kinematic hardening

In small strain theory, the use of Armstrong-Frederick-Chaboche models is widely accepted [29, 66, 67]. This model was enriched by a static recovery term, initially proposed by [30]. This term allows for a full or partial recovery of the kinematic hardening variable. It has been demonstrated that the Armstrong-Frederick model does not admit a dissipation potential [28] when used with a standard yield function [68]. Accordingly, this model is non standard. Generalized standard materials (GSM following [45]) are characterized by a single potential to describe the yield function, the flow rule and the evolution laws for internal variables. Nonlinear kinematic hardening can be introduced in this GSM framework by modifying the yield function as follows [68]

$$ f(\tilde{\mathbf{M}}, \tilde{\mathbf{X}}) = [\tilde{\mathbf{M}} - \tilde{\mathbf{X}}]_{eq} - R_0 + \frac{D}{2C} J^2(\tilde{\mathbf{X}}) \quad (28) $$

where $C$ and $D$ are material parameters related to kinematic hardening evolution and $J(\diamond) = \sqrt{\frac{3}{2}(\diamond : \diamond)}$. Nevertheless, this modification Eq. (28) induces an isotropic hardening term in addition to the one describing the kinematic hardening effect. A different approach has been proposed in [69] by introducing a new class of materials called implicit standard materials.

The extension of the Armstrong-Frederick model to finite strain ranges can be achieved in several ways. This issue was investigated thoroughly in [21, 32, 37]. A widely used approach is based on the multiplicative decomposition of the plastic part of the deformation gradient [7, 21, 32, 47]

$$ \tilde{\mathbf{F}}^p = \tilde{\mathbf{E}}_s^p \tilde{\mathbf{E}}_d^p \quad (29) $$

The free energy, depending on $\tilde{\mathbf{E}}^e$ and $\tilde{\mathbf{E}}_s^p = \frac{1}{2}(\tilde{\mathbf{F}}_s^{pT} \tilde{\mathbf{F}}_s^p - \tilde{\mathbf{1}})$, is split into elastic and kinematic hardening contributions as

$$ \rho_e \psi(\tilde{\mathbf{E}}^e, \tilde{\mathbf{E}}_s^p) = \rho_e \psi^e(\tilde{\mathbf{E}}^e) + \rho_e \psi^{\text{kin}}(\tilde{\mathbf{E}}_s^p) \quad (30) $$

The back stress acting on the local intermediate configuration is derived from

$$ \tilde{\mathbf{X}} = \rho_e \tilde{\mathbf{E}}_s^p \frac{\partial \psi^{\text{kin}}}{\partial \tilde{\mathbf{E}}_s^p} \tilde{\mathbf{E}}_s^{pT} \quad (31) $$

The dissipation inequality becomes

$$ \tilde{\mathbf{M}} : \tilde{\mathbf{D}}^p - \rho_e \frac{\partial \psi^{\text{kin}}}{\partial \tilde{\mathbf{C}}_s^p} : \dot{\tilde{\mathbf{C}}}_s^p \ge 0 \quad (32) $$

$$ \left( \tilde{\mathbf{M}} - 2\rho_e \tilde{\mathbf{E}}_s^p \frac{\partial \psi^{\text{kin}}}{\partial \tilde{\mathbf{C}}_s^p} \tilde{\mathbf{E}}_s^{pT} \right) : \tilde{\mathbf{D}}^p + 2\rho_e \tilde{\mathbf{C}}_s^p \frac{\partial \psi^{\text{kin}}}{\partial \tilde{\mathbf{C}}_s^p} : \tilde{\mathbf{D}}_d^p \ge 0 \quad (33) $$
where

$$ \mathcal{C}_s^p = \mathcal{F}_d^{p-T} \mathcal{C}_s^p \mathcal{F}_d^{p-1} \quad (34) $$

and

$$ \begin{aligned} \dot{\mathcal{C}}_s^p &= 2 \mathcal{F}_s^{pT} \bar{\mathcal{D}}^p \mathcal{F}_s^p - 2 \operatorname{sym}(\mathcal{C}_s^p \mathcal{D}_d^p), \\ \mathcal{D}_d^p &= \operatorname{sym}(\dot{\mathcal{C}}_d^p \mathcal{F}_d^{p-1}) \end{aligned} \quad (35) $$

We denote by $\tilde{\mathcal{X}} = 2\rho_e \mathcal{F}_s^p \frac{\partial \psi^{\text{kin}}}{\partial \mathcal{C}_s^p} \mathcal{F}_s^{pT}$ the back stress acting on the local intermediate configuration and by $\tilde{\mathcal{M}}_d = 2\rho_e \mathcal{C}_s^p \frac{\partial \psi^{\text{kin}}}{\partial \mathcal{C}_s^p}$ a Mandel-like stress tensor. The evolution equations satisfying (33) are given by

$$ \left\{ \begin{aligned} \dot{\mathcal{F}}_s^p \mathcal{F}_s^{p-1} &= \dot{\lambda} \frac{\partial f}{\partial \tilde{\mathcal{M}}}, \\ \dot{\mathcal{F}}_d^p \mathcal{F}_d^{p-1} &= \dot{\lambda} \frac{b}{c} \tilde{\mathcal{M}}_d \end{aligned} \right. \quad (36) $$

where b and c are material constants. For time-dependent plasticity, the plastic multiplier is given by $\dot{\lambda} = \frac{\partial \Omega}{\partial f}$ where $\Omega$ is the dissipation potential. Otherwise, the plastic multiplier $\dot{\lambda}$ is determined by applying the consistency condition $\dot{f} = 0$.

## The proposed kinematic hardening model

As investigated in [32, 37], different extensions of the Armstrong-Frederick model lead to qualitatively similar results if the material parameters are selected appropriately. Another formulation of the kinematic hardening model with both static and dynamic recovery terms is adopted here avoiding any further decomposition of $\mathcal{F}^p$ and resulting from the non-standard evolution equation as

$$ \tilde{\mathcal{X}} = \frac{2}{3} C \tilde{\alpha}, \quad \dot{\tilde{\alpha}} = \bar{\mathcal{D}}^p - \dot{\lambda} D \tilde{\alpha} - \frac{3}{2} \left( \frac{J(\tilde{\mathcal{X}})}{M} \right)^m \frac{\tilde{\mathcal{X}}}{J(\tilde{\mathcal{X}})} \quad (37) $$

where C, D, m and M are material parameters. The second term in the right describes the dynamic recovery and the third one is responsible for static recovery. The fact that the chosen internal variable $\tilde{\alpha}$ lives in the objective isoclinic intermediate configuration allows the use of the standard time derivative in the evolution equation. The push-forward operation applied to this equation will deliver the suitable objective derivative. This evolution equation must be complemented by the initial value of $\tilde{\alpha}$. The kinematic hardening variable $\tilde{\alpha}$ in Eq. (37) is symmetric provided that its initial value is symmetric.

This rule can be generalized to a non-symmetric kinematic hardening variable by substituting $\bar{\mathcal{D}}^p$ in Eq. (37) by $\bar{\mathcal{L}}^p$

$$ \dot{\tilde{\alpha}} = \bar{\mathcal{L}}^p - \dot{\lambda} D \tilde{\alpha} - \frac{3}{2} \left( \frac{J(\tilde{\mathcal{X}})}{M} \right)^m \frac{\tilde{\mathcal{X}}}{J(\tilde{\mathcal{X}})} \quad (38) $$

It follows that $\tilde{\mathcal{X}}$ is generally not symmetric, like the Mandel stress tensor. We mention also that a model with non-symmetric internal variable of kinematic hardening has been proposed in [46].

A more accurate description of a large variety of experimental stress-strain curves is possible by combining several independent kinematic hardening variables [30] such that

$$ \tilde{\mathcal{X}} = \sum_{i=0}^{N_x} \tilde{\mathcal{X}}_i, \quad \tilde{\mathcal{X}}_i = \frac{2}{3} C_{ij} \tilde{\alpha}_j \quad (39) $$

where $N_x$ is the total number of kinematic hardening variables associated with (visco)-plastic mechanisms. The matrix $C_{ij}$ accounts for the interaction between kinematic hardening variables.

## 2.6 Plastic spin

The tensor $\bar{\mathcal{L}}$, defined as the pull-back of $\bar{\mathcal{L}}$ to the intermediate configuration, can then be split into purely elastic and plastic parts as

$$ \bar{\mathcal{L}} = \mathcal{F}^{e-1} \bar{\mathcal{L}} \mathcal{F}^e = \bar{\mathcal{L}}^e + \bar{\mathcal{L}}^p \quad (40) $$

noting

$$ \bar{\mathcal{L}}^e = \mathcal{F}^{e-1} \dot{\mathcal{F}}^e \quad \text{and} \quad \bar{\mathcal{L}}^p = \dot{\mathcal{F}}^p \mathcal{F}^{p-1} \quad (41) $$

Strain rate and spin tensors expressed in the intermediate isoclinic configuration are derived as

$$ \bar{\mathcal{L}} = \bar{\mathcal{D}} + \bar{\mathcal{W}}, \quad \bar{\mathcal{L}}^e = \bar{\mathcal{D}}^e + \bar{\mathcal{W}}^e, \quad \bar{\mathcal{L}}^p = \bar{\mathcal{D}}^p + \bar{\mathcal{W}}^p \quad (42) $$

The plastic deformation is induced by various mechanisms such as mutual slips between material particles (crystals in metals, soil particles, etc.) without causing the rotation of the substructure defined by appropriate directors. Hence, the rotation and stretch of the substructure is induced only by the elastic distortion, including additional rigid body rotation, both included in $\mathcal{F}^e$. Accordingly, the spin of the substructure $\bar{\mathcal{W}}^e$ is independent of plastic deformation. Therefore, $\bar{\mathcal{W}}^e$ is given by subtracting the plastic spin from the total spin as follows

$$ \bar{\mathcal{W}}^e = \bar{\mathcal{W}} - \bar{\mathcal{W}}^p = \text{skw}(\mathcal{F}^{e-1} \dot{\mathcal{F}}^e) \quad (43) $$
One should note that if the plastic spin $\overline{\mathbf{W}}^p$ vanishes, the spin of directors will coincide with the material rotation rate.

Two important classes of plastic spin constitutive laws are distinguished in the literature. The first one is derived from normality rules whereas the second one follows from the application of tensor representation theorems.

### 2.6.1 Plastic spin derived from tensor representation theorems

A general explicit expression of plastic spin has been proposed independently by [40, 41] using appropriate tensor representation theorems as

$$ \overline{\mathbf{W}}^p = \beta_1(\underline{\mathbf{a}}\underline{\mathbf{s}} - \underline{\mathbf{s}}\underline{\mathbf{a}}) + \beta_2(\underline{\mathbf{a}}^2\underline{\mathbf{s}} - \underline{\mathbf{s}}^2\underline{\mathbf{a}}) + \beta_3(\underline{\mathbf{a}}\underline{\mathbf{s}}^2 - \underline{\mathbf{s}}^2\underline{\mathbf{a}}) \\ + \beta_4(\underline{\mathbf{a}}\underline{\mathbf{s}}\underline{\mathbf{a}}^2 - \underline{\mathbf{a}}^2\underline{\mathbf{s}}\underline{\mathbf{a}}) + \beta_5(\underline{\mathbf{s}}\underline{\mathbf{a}}\underline{\mathbf{s}}^2 - \underline{\mathbf{s}}^2\underline{\mathbf{a}}\underline{\mathbf{s}}) + \dots \quad (44) $$

where $\underline{\mathbf{s}}$ is the stress measure, $\underline{\mathbf{a}}$ denotes an internal variable and $\beta_i$ are material parameters. Tensor $\underline{\mathbf{a}}$ can be also regarded as a structure tensor [70, 71]. For instance, for a unidirectional composite described by a director $\underline{n}$ [72], $\underline{\mathbf{a}} = \underline{n} \otimes \underline{n}$. This allows to describe the evolution of the material substructure or directors [44]. To the best of our knowledge, all the studies of models including plastic spin make use only of the first order approximation of Eq. (44). By doing so, the plastic spin expressed in terms of Mandel's stress tensor is given by

$$ \overline{\mathbf{W}}^p = \beta(\underline{\mathbf{M}}\overline{\mathbf{D}}^p - \overline{\mathbf{D}}^p\underline{\mathbf{M}}) \quad (45) $$

This model is used to describe the evolution of anisotropic axes [42, 44] where $\beta$ is a parameter having the dimension of the inverse of stress. $\overline{\mathbf{W}}^p$ vanishes as long as $\underline{\mathbf{M}}$ and $\overline{\mathbf{D}}^p$ are coaxial. A similar but more sophisticated formulation was proposed in [73] to model deformation induced anisotropy in free-end torsion.

### 2.6.2 Plastic spin derived from the normality rule

According to this approach, the plastic spin is derived naturally from the normality rule as the skew-symmetric part of $\overline{\mathbf{L}}^p$

$$ \overline{\mathbf{W}}^p = \text{skw}(\overline{\mathbf{L}}^p) = \text{skw}\left(\frac{\partial\Omega(\underline{\mathbf{M}}, A_I)}{\partial\underline{\mathbf{M}}}\right) \quad (46) $$

This model of plastic spin will be studied in the case of crystal plasticity (see section 4.6).

## 3 Implementation in an object oriented FEM code

### 3.1 Object oriented architecture

The full description of the implementation of a generic formulation of constitutive equations in an object-oriented code, such as the Z-set software<fnref n="2" /> used in the present work, is detailed in [48, 68]. The software is designed in such a way that material models are implemented independently from the FEM. Therefore, the material library can be used by other FEM codes. This is made possible by proper interfaces between the FEM software and the material library Z-mat. The implementation of a generic material behaviour requires the following ingredients:

* **grad**: the imposed variable e.g. $\underline{\mathbf{F}}$. It allows driving the behaviour externally.

* **flux**: returned variable associated to a grad variable. It represents the response of the behavior law to the application of grad variable, e.g. the conjugate stresses $\underline{\mathbf{a}}$, $\underline{\mathbf{P}}$

* **EP**: external parameters as temperature, humidity, grain size, ... EP are set by the user and thus always known in advance.

* **IV**: Integrated variables which are to be integrated over a given time increment in order to update flux, e.g. $\underline{\mathbf{F}}^e$, $\underline{\mathbf{a}}_I$, ...

* **AV**: auxiliary variables do not define the material state directly, and they are kept for output. They may be useful in post-processing.

* **CO**: material parameters appearing as coefficients introduced in constitutive laws. They may depend on EP, IV and AV.

Figure 2 describes the implementation of a class named `GEN_FEFP` derived from `BEHAVIOR`. This class allows to consider various inelastic mechanisms through the class `POTENTIAL`. Each potential includes a flow rule (time-independent plasticity, Norton power law...), several isotropic and kinematic hardening rules. The internal variables associated with kinematic hardening (i.e. $\underline{\mathbf{a}}$) are held by the corresponding `KINEMATIC` object. Explicit and implicit integration of constitutive equations are handled by `RUNGE_KUTTA` and `THETA_METHOD` classes, respectively. The `CRITERION` object specifies the yield surface (von Mises, Tresca, Hill,...). This implementation allows a minimum programming effort since all the required classes to build a material behavior (except `GEN_FEFP`) are already implemented and used in the framework of other formulations (small strain, hypoelasto-plasticity,...).

<footnote n="2">² See www.zset-software.com.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2 A diagram showing the organization of some objects used to build material behaviors</figcaption></figure>

## 3.2 Global resolution of equilibrium equations

The current (resp. reference) configuration of the body at time $t$ (resp. $t_0$) is called $V$ (resp. $V_0$) with boundary $\partial V$ (resp. $\partial V_0$). The latter can be split into the sub-boundaries $\partial V^u$ and $\partial V^{tr}$ such that $\partial V = \partial V^u \cup \partial V^{tr}$ and $\partial V^u \cap \partial V^{tr} = \emptyset$, where Dirichlet and Neumann conditions are respectively prescribed. Corresponding surfaces $\partial V_0^u$ and $\partial V_0^{tr}$ are defined on the reference boundary of the body. The space of kinematically admissible displacements field is the set of sufficiently regular displacement functions that satisfy the Dirichlet conditions:

$$ \mathcal{K} = \{ \underline{\mathbf{u}}(\underline{\mathbf{x}}) | \underline{\mathbf{u}} = \tilde{\underline{\mathbf{u}}}(\underline{\mathbf{x}}) \text{ if } \underline{\mathbf{x}} \in \partial V^u \} \quad (47) $$

knowing that $\tilde{\underline{\mathbf{u}}}$ is the prescribed displacement field on $\partial V^u$. Virtual displacements are set to zero over $\partial V^u$ defining the space of virtual displacements as

$$ \mathcal{V} = \{ \underline{\eta}(\underline{\mathbf{x}}) | \underline{\eta}(\underline{\mathbf{x}}) = 0 \text{ if } \underline{\mathbf{x}} \in \partial V^u \} \quad (48) $$

The surface tractions are prescribed over the region $\partial V^{tr}$.

Following [74] among others, the initial boundary value problem amounts to find a displacement field $\underline{\mathbf{u}} \in \mathcal{K}$ that satisfies

$$ \mathcal{W}(\underline{\mathbf{u}}, \underline{\eta}) = 0, \quad \forall \underline{\eta} \in \mathcal{V} \quad (49) $$

where the virtual work functional under finite deformation is given by

$$ \mathcal{W}(\underline{\mathbf{u}}, \underline{\eta}) = \int_{V_0} (\underline{\mathbf{P}} : \nabla_X \underline{\eta} - \underline{b} \cdot \underline{\eta}) dV - \int_{\partial V_0^{tr}} \underline{t} \cdot \underline{\eta} dS \quad (50) $$

where $\underline{b}$ and $\underline{t}$ denote the reference body force and surface traction fields, respectively. $\underline{\mathbf{P}}$ stands for the Boussinesq stress tensor. The Boussinesq stress tensor $\underline{\mathbf{P}}$ is related to Kirchhoff's stress tensor $\underline{\tau}$ by

$$ \underline{\mathbf{P}} = \underline{\tau} \underline{\mathbf{F}}^{-T} = J \underline{\sigma} \underline{\mathbf{F}}^{-T} \quad (51) $$

The linearization of Eq. (49) at a given state defined by the field $\underline{\mathbf{u}}^*$ is written

$$ \mathcal{W}(\underline{\mathbf{u}}^*, \underline{\eta}) + D\mathcal{W}(\underline{\mathbf{u}}^*, \underline{\eta})[\Delta\underline{\mathbf{u}}] = 0, \quad \forall \underline{\eta} \in \mathcal{V} \quad (52) $$

where $D\mathcal{W}(\underline{\mathbf{u}}^*, \underline{\eta})[\Delta\underline{\mathbf{u}}]$ is the directional derivative of $\mathcal{W}(\underline{\mathbf{u}}^*, \underline{\eta})$ in the direction of $\Delta\underline{\mathbf{u}}$. For convenience, the force and surface traction fields on $\partial V_0^{tr}$ are assumed to be independent of the displacement field. Then,

$$ D\mathcal{W}(\underline{\mathbf{u}}^*, \underline{\eta})[\Delta\underline{\mathbf{u}}] = \left. \frac{d}{d\epsilon} \left( \mathcal{W}(\underline{\mathbf{u}}^* + \epsilon \Delta\underline{\mathbf{u}}, \underline{\eta}) \right) \right|_{\epsilon=0} \quad (53) $$

$$ = \frac{d}{d\epsilon} \left( \int_{V_0} \underline{\mathbf{P}}(\underline{\mathbf{F}}^* + \epsilon \Delta\underline{\mathbf{F}}) : \nabla_X \underline{\eta} \right) \Big|_{\epsilon=0} \quad (54) $$

$$ = \frac{d}{d\epsilon} \left( \int_{V_0} [\underline{\mathbf{P}}(\underline{\mathbf{F}}^*) + \Delta\underline{\mathbf{P}}(\underline{\mathbf{F}}^*, \epsilon \Delta\underline{\mathbf{F}})] : \nabla_X \underline{\eta} \right) \Big|_{\epsilon=0} \quad (55) $$

$$ = \int_{V_0} \left[ \left. \frac{\partial \underline{\mathbf{P}}}{\partial \Delta \underline{\mathbf{F}}} \right|_{\underline{\mathbf{F}}^*} : \Delta \underline{\mathbf{F}} \right] : \nabla_X \underline{\eta} \quad (56) $$

where

$$ \Delta \underline{\mathbf{F}} = \nabla_X (\Delta \underline{\mathbf{u}}), \quad \underline{\mathbf{F}}^* = \underline{\mathbf{1}} + \nabla_X \underline{\mathbf{u}}^* \quad (57) $$

The tangent modulus is computed as

$$ \underline{\mathcal{A}} = \left. \frac{\partial \underline{\mathbf{P}}}{\partial \Delta \underline{\mathbf{F}}} \right|_{\underline{\mathbf{F}}^*} = \left. \frac{\partial \underline{\tau}}{\partial \Delta \underline{\mathbf{F}}} \right|_{\underline{\mathbf{F}}^*} \underline{\mathbf{F}}^{-T} + \underline{\tau} \left. \frac{\partial \underline{\mathbf{F}}^{-T}}{\partial \Delta \underline{\mathbf{F}}} \right|_{\underline{\mathbf{F}}^*} \quad (58) $$

or in index notation

$$ \mathcal{A}_{ijkl} = \frac{\partial \tau_{ip}}{\partial \Delta F_{kl}} F_{jp}^{-1} - \tau_{ip} F_{jk}^{-1} F_{lp}^{-1} \quad (59) $$

By approximating the domain $V_0$ with a finite number of elements $n_{el}$ denoted by $V_0^e$, the discrete form of the virtual work (50) is written

$$ \sum_{e=0}^{n_{el}} \int_{V_0^e} ( [\underline{\mathbf{B}}]^T \{\underline{\mathbf{P}}\} - [\underline{\mathbf{N}}^T] \{\underline{\mathbf{b}}\} ) dV $$
$$-\sum_{e=0}^{n_{el}} \int_{\partial V_0^{tr,e}} [\mathbf{N}]^T \{\mathbf{t}\} dS = 0 \qquad (60)$$

$[\mathbf{N}]$ and $[\mathbf{B}]$ denote the interpolation matrix and the discrete material gradient operator, respectively (see Appendix B). The element stiffness matrix is given by

$$[\mathbf{K}^e] = \int_{V_0^e} [\mathbf{B}]^T [\mathbf{A}] [\mathbf{B}] dV \qquad (61)$$

The global stiffness matrix is obtained by assembling the element tangent stiffness matrices as

$$[\mathbf{K}^g] = \mathbf{A}_{e=1}^{n_{el}} ([\mathbf{K}^e]) \qquad (62)$$

The linear set of equations to be solved for $\Delta\underline{\mathbf{u}}$ iteratively is given by

$$[\mathbf{K}]^g \Delta\underline{\mathbf{u}} = -\mathcal{R}(\underline{\mathbf{u}}) \qquad (63)$$

where

$$\mathcal{R}(\underline{\mathbf{u}}) = f^{\text{int}}(\underline{\mathbf{u}}) - f^{\text{ext}} \qquad (64)$$

The global internal and external forces are written

$$f^{\text{int}} = \mathbf{A}_{e=1}^{n_{el}} \left( \int_{V_0^e} [\mathbf{B}]^T \{\mathbf{P}\} dV \right) \qquad (65)$$

$$f^{\text{ext}} = \mathbf{A}_{e=1}^{n_{el}} \left( \int_{V_0^e} [N^T]\{\mathbf{b}\} dV + \int_{\partial V_0^{tr,e}} [N]^T \{\mathbf{t}\} dS \right) \qquad (66)$$

The consistent tangent moduli $\mathcal{A}^c$ are calculated from the implicit incremental constitutive equations as

$$\frac{\partial \underline{\mathbf{\tau}}(\alpha_n, \mathbf{F}_{n+1})}{\partial \Delta \mathbf{F}_{n+1}} = \frac{\partial \underline{\mathbf{\tau}}(\alpha_n, \mathbf{F}_{n+1})}{\partial \Delta \mathbf{F}_{n+1}^e} \frac{\partial \Delta \mathbf{F}_{n+1}^e}{\partial \Delta \mathbf{F}_{n+1}} \qquad (67)$$

The general iterative resolution algorithm of the incremental boundary value problem is summarized in Fig. 3.

## 3.3 Integration of constitutive equations

The set of time-integrated variables is given by

$$\mathcal{V}_{\text{int}} = \{ \mathcal{F}^e; p_1, \alpha_1; \dots; p_n, \alpha_n \} \qquad (68)$$

$p_i$ is the accumulated inelastic strain associated with the ith mechanism, $\alpha_i$ denote internal variables describing both isotropic (e.g. $r_i$) and kinematic hardening (e.g. $\alpha_i$) specific to each mechanism. Semicolons in (68) represent the separation between different mechanisms. A distinction between $\alpha_i$ and

**Result:** Compute the displacement $\underline{\mathbf{u}}_{n+1}$ at $t_{n+1}$

(i) Set initial guess $\underline{\mathbf{u}}_{n+1}^{(0)}$, $k=1$

while $k < N_{\text{iter}}^{\text{max}}$ do

(ii) $\underline{\mathbf{u}}_{n+1}^{(k)} = \underline{\mathbf{u}}_{n+1}^{(k-1)} + \Delta\underline{\mathbf{u}}^{(k)}$

(iii) $\underline{\mathbf{F}}_{n+1}^{(k)} = \underline{\mathbf{z}} + \nabla \chi \underline{\mathbf{u}}_{n+1}^{(k)}$

(iv) Update stress and state variables using constitutive equations

(v) Compute consistent tangent moduli $\mathcal{A}^c$

(vi) Calculate residuals $\mathcal{R}(\underline{\mathbf{u}}_{n+1}^{(k)}) = f^{\text{int}}(\underline{\mathbf{u}}_{n+1}^{(k)}) - f_{n+1}^{\text{ext}}$

if convergence then

go to next increment $n+1$ (i)

else

(vii) Compute global stiffness matrix $[\mathbf{K}]^g$

(viii) $\Delta\underline{\mathbf{u}}^{(k+1)} = -[\mathbf{K}]^g^{-1}\mathcal{R}$, $k \leftarrow k+1$

end

end

<caption>**Fig. 3** Iterative resolution algorithm of the incremental boundary value problem</caption>

$p_i$ was considered so that partial derivatives required by integration methods can be implemented efficiently. The model for each deformation mechanism will be defined by the following system of elementary equations via the POTENTIAL interface:

* Plasticity criteria $f_i(\underline{\mathbf{M}}, p_i, \alpha_i)$: Each criterion depends on the Mandel stress tensor, the accumulated plastic strain and variables that describe multi-kinematic hardening. Only the case where a criterion is associated with one and only one dissipation potential is discussed in the following. Alternative models involving multi-mechanisms and one single unified criterion can be found in [65]. The implementation of this kind of models is discussed thoroughly in [48].

* Flow rule: The plastic multiplier is defined in two ways. For time-independent plasticity, $\dot{\lambda}$ must fulfill the consistency condition Eq. (17). In the rate-dependent case, the multiplier is defined as the derivative of dissipation potential $\Omega_i$ w.r.t the yield function $f_i$ Eq. (18). It follows that

$$\overline{\mathbf{N}}_i^p = \dot{\lambda}_i \frac{\partial f_i}{\partial \overline{\mathbf{N}}} = \dot{\lambda}_i \mathcal{N}_i \qquad (69)$$

with $\mathcal{N}_i$ being the inelastic flow direction.

* Isotropic and kinematic hardening laws take the generic form

$$\dot{\alpha}_i = \dot{\lambda}_i \mathbf{m}_i - \dot{\mathbf{q}}_i \qquad (70)$$

where $\mathbf{m}_i$ is the hardening potential normal and $\dot{\mathbf{q}}_i$ represents the time derivative of the hardening variable evolution due to static recovery effects. Note that in case of associated plasticity, $\mathbf{m}_i = \frac{\partial f_i}{\partial \mathbf{A}_i}$, meaning that the hardening potential normal is simply given by the normal to the yield surface. For instance, $\underline{\mathbf{m}}_i = \underline{\mathbf{N}}_i - \frac{3D_i}{2C_i} \underline{\mathbf{X}}_i$
for nonlinear kinematic hardening. The evolution of the nonlinear isotropic hardening variable $r$, from Eq. (14), writes $\dot{r} = \dot{\lambda}(1 - \frac{R}{Q})$ which implies that $m = 1 - \frac{R}{Q}$. The integrated variable $r$ is related to $R$ by $R = b\dot{Q}r$ (see Eq. (15)), with the material parameters $b$ and $Q$. The latter model can be integrated analytically (cf. Equation (104) in Sect. 5), which allows to reduce the number of time-integrated variables in the code.

Two integration methods for ordinary differential equations have been implemented for the present general formulation. The first one is explicit, namely second or fourth order Runge-Kutta methods with automatic time-stepping, the second one is implicit: $\theta$-method resolved using the iterative Newton-Raphson scheme.

Each inelastic deformation mechanism is accounted for within the class POTENTIAL. This class provides the increment of integrated variables for explicit integration and residuals for implicit integration. Each POTENTIAL object possesses various methods with regard to the model definition (isotropic and kinematic hardening, yield criteria,...). The implementation supports an unlimited number of POTENTIAL objects with possible interactions. In case of interactions, another class is dedicated to add interaction terms appropriately. The expression of $\tilde{\mathcal{L}}^p$ is calculated through summation of all potential contributions according to Eq. (27), in addition to the supplementary ad-hoc constitutive equation for the plastic spin.

**RUNGE_KUTTA:**

The integration of variables is based on the calculation of the rate of $V_{int}$. For viscoplastic cases, the plastic multiplier is calculated using the flow rule. For time-independent plasticity, the increment of the plastic multiplier is derived from the consistency condition as

$$2\tilde{N} : \left[ \mathbf{1} \otimes \mathbf{\Pi}^e + \frac{1}{2} \tilde{\mathbf{C}}^e \cdot \mathbf{\mathcal{C}} \right] : \left( \tilde{\mathbf{F}}^{eT} \tilde{\mathbf{D}} \tilde{\mathbf{F}}^e - \tilde{\mathbf{C}}^e \sum_{i}^{N} \tilde{\mathcal{L}}_{\tilde{i}}^p \right) + \frac{\partial f}{\partial \tilde{A}_i} : \sum_{i}^{N} \frac{\partial \tilde{A}_i}{\partial \alpha_i} \dot{\alpha}_i + \frac{\partial f}{\partial p_i} \dot{p}_i + \frac{\partial f}{\partial \text{EP}} : \dot{\text{EP}} = 0 \quad (71)$$

**THETA_METHOD:**

$\theta = 0$ corresponds to the explicit Euler method and $\theta = 1$ results in the so-called backward Euler method. The constitutive equations are expressed in the residual form as follows

$$\mathcal{R}_{el} = \tilde{\mathbf{F}}^e \tilde{\mathbf{F}}^p - \tilde{\mathbf{F}} \quad \text{or} \quad \mathcal{R}_{el} = \Delta \tilde{\mathbf{F}}^e - \Delta \tilde{\mathbf{F}} \tilde{\mathbf{F}}^{-1} \tilde{\mathbf{F}}_{cor}^e + \tilde{\mathbf{F}}_{cor}^e \sum_{i=0}^{N} \Delta p_i \tilde{N}_i \quad (72)$$

$\mathcal{R}_{p_i} = f_i(\tilde{M}, \tilde{A}_i)$ or

$$\mathcal{R}_{p_i} = \Delta p_i - \Delta t \frac{\partial \Omega}{\partial f_i} \quad (\text{viscoplasticity}) \quad (73)$$

$$\mathcal{R}_{\alpha_i} = \Delta \alpha_i - \Delta p_i \mathbf{m}_i - \Delta t \dot{\mathbf{q}}_i \quad (74)$$

$N$ being the total number of mechanisms. The residual in Eq. (72)-left involves the use of the exponential mapping which satisfies the plastic incompressibility [75, 76]. Its linearization requires the expression of the derivative of the exponential of a second order tensor w.r.t. a second order tensor. As $\tilde{\mathcal{L}}^p$ is in general non symmetric, calculating these terms is a non-trivial task. In that case, the infinite series representation is used and truncated [32, 47]. Alternatively, in Eq. (72)b, the plastic incompressibility is not satisfied anymore. For this reason, the elastic part of the deformation gradient is corrected at each iteration as

$$\tilde{\mathbf{F}}_{cor}^e = \left( \frac{\det \tilde{\mathbf{F}}}{\det \tilde{\mathbf{F}}^e} \right)^{\frac{1}{3}} \tilde{\mathbf{F}}^e \quad (75)$$

The values of all associated forces and parameters, calculated from internal variables evaluated at an intermediate time designated by $\theta$ are

$$\mathcal{V}_{int}^{t+\theta\Delta t} = \mathcal{V}_{int}^{t} + \theta \Delta \mathcal{V}_{int} \quad (76)$$

The set of equations (72,73,74) can be gathered in the following form

$$\mathcal{R}(\mathcal{V}_{int}^{t+\theta\Delta t}, \Delta\mathcal{V}_{int}) = 0 \quad (77)$$

Since Eq. (77) is highly non-linear, it is usually solved by means of a Newton method which requires the calculation of the Jacobian matrix [J]. The new estimate of $\Delta\mathcal{V}_{int}^{k+1}$ is then given by

$$\Delta\mathcal{V}_{int}^{k+1} = \Delta\mathcal{V}_{int}^{k} - [J]^{-1}\mathcal{R}^{k} \quad (78)$$

where

$$[J] = \left. \frac{\partial \mathcal{R}}{\partial \Delta \mathcal{V}_{int}} \right|_{k+1} \quad (79)$$

and $\mathcal{R}^k$ denotes the local residual at the $k$-th iteration. The variation of $\mathcal{R}$ resulting from the variation of $\mathcal{V}_{int}$ and $\tilde{\mathbf{F}}$ vanishes as well

$$\delta \mathcal{R} = \frac{\partial \mathcal{R}}{\partial \Delta \mathcal{V}_{int}} \delta \Delta \mathcal{V}_{int} + \frac{\partial \mathcal{R}}{\partial \Delta \tilde{\mathbf{F}}} \delta \Delta \tilde{\mathbf{F}} = 0 \quad (80)$$

which implies that, after convergence, the inverted Jacobian matrix relates the change of $\Delta\mathcal{V}_{int}$ with respect to a change in $\Delta\tilde{\mathbf{F}}$ as
<figure><img src="image_4.png" /><figcaption>**Fig. 4** Diagram showing explicit and implicit integration methods. Each POTENTIAL provides the update of associated internal variables and the corresponding block of the Jacobian matrix</figcaption></figure>

$$ \delta\Delta V_{int} = - \left( \frac{\partial \mathcal{R}}{\partial \Delta V_{int}} \right)^{-1} \frac{\partial \mathcal{R}}{\partial \Delta F} \delta\Delta \mathcal{F} \quad (81) $$

The expression of the Jacobian matrix is detailed in the Appendix (C). By using Eq. (81), one can calculate the term $\frac{\partial \Delta \mathcal{F}^e}{\partial \Delta \mathcal{F}}$ in Eq. (67).

## 4 Applications to volume element simulations

The capabilities of the model formulation and its implementation are illustrated in this section in the case of complex homogeneous loading conditions. The tests are therefore performed at the material point level and require integration of the constitutive equations. Some original features of the models are highlighted and compared to the predictions of hypoelastic models.

### 4.1 Cyclic closed deformation path

Hypoelastic formulations are well known to result in some spurious predictions under complex loading conditions. For instance, when a closed strain cycle is applied, the resulting stress cycle is not necessarily closed. This issue was investi-

gated theoretically [5, 6] and numerically [7] in the case of purely mechanical loadings. In this example, the response of finite strain formulations is investigated in the thermoelastic domain for two different loading cases. In both cases, we consider two hypoelastic constitutive models for the Cauchy stress based on the Jaumann and Green–Naghdi rates and a hyperelastic model given by Eq. (23).

**Case 1:**

A cyclic and non-proportional mechanical loading is applied (see Fig. 5a). Accordingly, the applied deformation gradient and the resulting Cauchy stress tensor have the following forms

$$ \begin{pmatrix} F_{11}(t) & F_{12}(t) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \quad \text{and} \quad \boldsymbol{\sigma} = \begin{pmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_{22} & 0 \\ 0 & 0 & \sigma_{33} \end{pmatrix} \quad (82) $$

where the functions $F_{11}(t)$ and $F_{12}(t)$ are prescribed.

**Case 2:**

Is concerned with a non-proportional thermomechanical loading described in Fig. 5b and corresponding to prescribed shear $F_{12}(t)$ and temperature variation $\Delta T(t)$. The deformation gradient is imposed as
<figure><img src="image_5.png" /><figcaption>**Fig. 6** Case 1, the resulting stress during the first and the 10000th cycles for **a** hyperelastic material **b** hypoelastic materials. Young's modulus $E = 210000$ MPa, and Poisson ratio $\nu = 0.3$</figcaption></figure>

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

$$ \tilde{F} = \begin{pmatrix} 1 & F_{12}(t) & 0 \\ 0 & F_{22}^* & 0 \\ 0 & 0 & F_{33}^* \end{pmatrix} \quad \text{and} \quad \tilde{\sigma} = \begin{pmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \quad (83) $$

$F_{22}^*$ and $F_{33}^*$ are determined by the analysis in such a way that the conjugate Cauchy stress components vanish. The component $\sigma_{11}$ does not vanish due to the applied temperature and to the Poynting effect, i.e. finite deformation induced stress components other than shear. The Figs. 6 and 7 depict an elastically inconsistent response of the two hypoelastic formulations, namely Jaumann and Green–Naghdi. When a closed strain cycle is applied, residual stresses remain at the end of each cycle. Consequently, due to the accumulation of residual stresses during the deformation process, the resulting stress drifts away over cycles. In contrast, for a hyperelastic model, no residual stresses are detected after each cycle, i.e. all stress components return back to zero. As suggested by [8], some hypoelastic models, e.g. based on the logarithmic rate, produce consistent results compared to hyperelastic models. Nevertheless, in the case of elastoplasticity, these rates are no longer integrable [10]. It has been shown that any hypoelastic law is integrable in the case of proportional loading i.e. for deformation processes depending on a single parameter [77].

## 4.2 Simple glide with kinematic hardening

The deformation gradient for a simple glide test has the form

$$ \tilde{F} = \begin{pmatrix} 1 & \gamma(t) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \quad (84) $$

In this section, a comparison is drawn between the elastoplastic models previously presented by looking at the classical cases of monotonic and cyclic simple shear. In particular the proposed kinematic hardening evolution law (37) is compared to the existing formulation (31), for which the kinematic hardening contribution $\psi^{\text{kin}}$ to the free energy is given in the Neo-Hookean form [32, 47]

$$ \psi^{\text{kin}}(\tilde{E}_s^p) = c \left( \text{tr}(\tilde{E}_s^p) - \log(\det \tilde{F}_s^p) \right) \quad (85) $$

In the case of linear kinematic hardening, significant differences are found between the various formulations as shown in Figs. 8 and Fig. 9. Figure 8 depicts a spurious oscillatory response of both multiplicative and Jaumann formulations. The corresponding analytical solution in the rigid plastic case can be found in [68]. In contrast, the model based on the Green–Naghdi rate predicts that $\sigma_{12}$ increases monotonically. The model based on the multiplicative decomposition of the plastic part of the deformation gradient does not display oscillations as depicted in Fig. 9a. These oscillations
<figure><img src="image_7.png" /><figcaption>**Fig. 8** Monotonic simple glide test for different formulations of elastoplasticity with kinematic hardening: Stress oscillations in the case of Jaumann rate and the multiplicative decomposition if $D = 0$. These oscillations are suppressed as the parameter $D$ increases. Values of material parameters: $C = 10000$ MPa, $R_0 = 1000$ MPa</figcaption></figure>

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

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

can be suppressed by increasing the value of the parameter $D$ from Eq. (37). By doing so, kinematic hardening saturates rapidly which leads to almost the same predictions by different models. Accordingly, when the dynamic recovery term is sufficiently high (compared to the storage part), the saturation rate of the various models becomes similar.

### 4.3 Static recovery of kinematic hardening

In this example, the effect of the static recovery of the kinematic hardening is illustrated. It corresponds to the last term in the evolution equation Eq. (37). This term introduces time-dependent material behavior even in the absence of viscosity. We consider a von Mises surface yield given by

$$f(\underline{\underline{M}}, R) = ((\underline{\underline{M}} - \underline{\underline{X}})^{dev} : (\underline{\underline{M}} - \underline{\underline{X}})^{dev})^{1/2} - R_0 \quad (86)$$

where $(\diamond)^{dev}$ denotes the deviatoric part. The back stress $\underline{\underline{X}} = \frac{2}{3}C\underline{\underline{\alpha}}$ and the evolution of $\underline{\underline{\alpha}}$ follows the constitutive equation (37). This model is applied to relaxation and creep tests.

#### 4.3.1 Relaxation test under simple glide

A time-dependent simple glide function $\gamma(t)$ is considered according to the following conditions

$$\underline{\underline{F}} = \begin{pmatrix} 1 & \gamma(t) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \gamma(t) = \begin{cases} \frac{t}{2t_0} & 0 \le t \le 1 \text{ s} \\ \frac{1}{2} & 1 \text{ s} \le t \le 4 \text{ s} \end{cases} \quad (87)$$

with $t_0 = 1$ s. This loading corresponds to monotonic glide followed by a constant shear value inducing stress relaxation in order to highlight the impact of the static recovery term. The material response is shown in Fig. 10 for three different sets of values of material parameters $(m, M)$.

The following features can be observed

* $0 \le t \le 1$ s: During this time interval, the applied shear is monotonically increasing. After yielding, there is a competition between different terms in Eq. (37): storage part, dynamic and static recoveries. For given parameters $C$ and $D$, a higher value of parameter $M$ or a lower value of parameter $m$ allow for slower recovery of the kinematic hardening leading to a higher value of $\sigma_{12}$ (see Fig. 10).

* $1 \le t \le 4$ s: The imposed shear strain is constant i.e. $\dot{\underline{\underline{F}}} = 0$. Consequently, the variation of the accumu-
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Cyclic simple glide loading: **a** linear (D = 0) **b** nonlinear kinematic hardening with dynamic recovery</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Relaxation test: Higher values of *m* and lower values of *M* result in a faster recovery of kinematic hardening. *C* = 300000 MPa, *D* = 20. Parameter *M* in MPa</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Plastic flow during a creep test</figcaption></figure>

lated plastic strain is low due to small elastic strain since $\bar{\underline{\zeta}}^p = -\underline{\zeta}^{e-1} \dot{\underline{\zeta}}^e$ during a relaxation test. Therefore, the evolution of the kinematic hardening variable reduces to

$$ \dot{\underline{\alpha}} \approx -\frac{3}{2} \left( \frac{\underline{J}(\underline{\underline{X}})}{M} \right)^m \frac{\underline{\underline{X}}}{\underline{J}(\underline{\underline{X}})} \quad (88) $$

The recovery rate increases then with higher values of the power *m* (resp. lower values of *M*).

In metals, this effect occurs significantly at high temperature due to thermal activation. In fact, the crystalline structure of the metal is partially recovered by annihilation of dislocations and redistribution of point defects [30]. This relaxation of internal stresses generally results in a decrease of the mechanical resistance.

### 4.3.2 Creep test

In this example, the imposed deformation components are $F_{ij} = 0$ for $i \neq j$ and $\{ij\} \neq \{12\}$. $F_{11}$, $F_{22}$, $F_{33}$ and $F_{12}$ are set to be free. The imposed Cauchy stress component $\sigma_{12}$ is

is a function of time (s) given by

$$ \sigma_{12} = \begin{cases} 1000\bar{t} \text{ MPa}, & 0 \le \bar{t} \le 1 \\ 1000 \text{ MPa}, & 1 \le \bar{t} \le 4 \end{cases} \quad (89) $$

where $\bar{t} = t/t_0$, $t_0 = 1$ s. The remaining components of the Cauchy stress tensor vanish.

The material response is shown in Fig. 11. It is important to note that the used model does not include any viscosity, meaning that the yield condition is exactly fulfilled. The creep deformation is induced solely by the static relaxation of the back-stress that has developed during the loading stage (see Eq. (37)):

* $0 \le \bar{t} \le 1$: As the imposed stress increases, the accumulated plastic strain increases after yielding ($\sigma_{12} = \sqrt{3}R_0$).

* $1 \le \bar{t} \le 4$: According to the consistency condition, the kinematic hardening does not evolve since the applied stress is constant. Therefore, the plastic multiplier is constant and the accumulated plastic multiplier increases linearly. Due the static recovery term, the model becomes time-dependent. In other words, if the static recovery term is omitted, then $\dot{\underline{\alpha}} = 0 \Rightarrow \dot{p} = 0$.
<table><caption>**Table 1** Material parameters for a model with two inelastic mechanisms: (p) time-independent plastic and (v) viscoplastic. $C_{pv} = C_{vp} = 0$ MPa</caption><thead><tr><th>E<br>(MPa)</th><th>ν<br>(-)</th><th>K<br>(MPa<sup>1/2</sup>)</th><th>n<br>(-)</th><th>R<sub>0v</sub><br>(MPa)</th><th>C<sub>v</sub><br>(MPa)</th><th>D<sub>v</sub><br>(-)</th><th>R<sub>0p</sub><br>(MPa)</th><th>C<sub>p</sub><br>(MPa)</th><th>D<sub>p</sub><br>(-)</th></tr></thead><tbody><tr><td>210000</td><td>0.3</td><td>120</td><td>7</td><td>0</td><td>20000</td><td>200</td><td>140</td><td>1000</td><td>10</td></tr></tbody></table>

## 4.4 Application to a von Mises-based multimechanism model

The concept of multimechanism modeling is applied in this section to isotropic von Mises plasticity. A model is introduced involving two inelastic mechanisms and two plasticity criteria, called 2M2C in the terminology defined in the references [64, 65]. The yield function, flow rule and evolution equations adopted for this example are as follows

$$f_p(\underline{M}, \underline{X}_p) = ((\underline{M} - \underline{X}_p)^{\text{dev}} : (\underline{M} - \underline{X}_p)^{\text{dev}})^{1/2} - R_{0p} \quad (90)$$

$$f_v(\underline{M}, \underline{X}_v) = ((\underline{M} - \underline{X}_v)^{\text{dev}} : (\underline{M} - \underline{X}_v)^{\text{dev}})^{1/2} - R_{0v} \quad (91)$$

$$\dot{v} = \left( \frac{f_v}{K} \right)^n, \\ \dot{p} = \frac{\underline{N}_p : (\underline{1} \otimes \underline{\Pi}^e + \frac{1}{2} \underline{C}^e \cdot \underline{\mathbb{C}}) : (\underline{F}^e T \underline{D} \underline{F}^e - \dot{v} \underline{C}^e \underline{N}_v)}{\underline{N}_p : \left[ (\underline{1} \otimes \underline{\Pi}^e + \frac{1}{2} \underline{C}^e \cdot \underline{\mathbb{C}}) : (\underline{C}^e \underline{N}_p) \right] + \frac{1}{3} C_p \underline{N}_p : \underline{m}_p} \quad (92)$$

$$\dot{\underline{\alpha}}_p = \dot{p} \left( \frac{\partial f_p}{\partial \underline{M}} - D_p \underline{\alpha}_p \right) \cdot \dot{\underline{\alpha}}_v = \dot{v} \left( \frac{\partial f_v}{\partial \underline{M}} - D_v \underline{\alpha}_v \right) \quad (93)$$

$$\left(\frac{\underline{X}_p}{\underline{X}_v}\right) = \frac{2}{3} \begin{pmatrix} C_p & C_{pv} \\ C_{vp} & C_v \end{pmatrix} \begin{pmatrix} \underline{\alpha}_p \\ \underline{\alpha}_v \end{pmatrix} \quad (94)$$

The first plastic mechanism associated with the yield function $f_p$ is rate-independent, whereas the second one is viscoplastic. The (visco)plastic multipliers are computed either by a power law or the consistency condition, according to Eq. (92). The material parameters used in the following examples are given in Table 1.

Consider now a simple tensile/compressive test under strain control, divided into five stages. The resulting stress state is uni-axial ($\sigma_{11}$). The response of the model for two distinct strain rates is given in Fig. 12. The following observations can be made:

* $F_{11} = 1.1t/t_0$, for $0 \le t \le t_1$ (tension): The (v) mechanism is active first because the corresponding yield stress is taken as $R_{0v} = 0$. But after a while, the mechanism (p) is activated once the threshold $R_p$ is reached. In fact, the activation of an inelastic mechanism will depend on the associated initial yield stress and also on the hardening properties of the other mechanism. Higher strain rates

induce more plastic strain, while lower strain rates result in more viscoplastic strain.

* $F_{11} = 1.1$, for $t_1 \le t \le t_2$: During this stage the imposed strain is maintained at a constant value. The inelastic deformation remains quasi-constant accompanied by stress relaxation due to dynamic recovery of kinematic hardening.

* $F_{11} = 1.1 - 0.1(t - t_2)/t_0$, for $t_2 \le t \le t_3$ (unloading+compression): The plastic yielding in compression occurs at a rather large stress level due to the Bauschinger effect.

* $F_{11} = 0.9$, for $t_3 \le t \le t_4$: The imposed strain is constant. A quasi-constant inelastic deformation and stress relaxation due to kinematic hardening recovery are observed.

* $F_{11} = 0.9 + 0.1(t - t_4)/t_0$, for $t_4 < t < t_5$ (unloading phase). At $F_{11} = 1$, the residual stress does not vanish due to kinematic hardening.

## 4.5 Plastic spin in anisotropic plasticity

The effect of plastic spin on the response of a Hill perfectly plastic material is illustrated in this example. Hill's yield criterion is expressed in terms of Mandel stress tensor in the form

$$f(\underline{M}) = \left( F(M_{22} - M_{33})^2 + G(M_{33} - M_{11})^2 + H(M_{11} - M_{22})^2 + 2NM_{12}^2 + 2LM_{23}^2 + 2MM_{13}^2 \right)^{1/2} - R_0 \quad (95)$$

where $F, G, H, N, M, L$ are material parameters characterizing the anisotropy of plasticity. This yield function reduces to von Mises if $F = G = H = 0.5$ and $H = N = M = 1.5$.

A simple glide test is considered with constant shear rate in the plane (1, 2). When the hypoelastic-plastic model is used with the Jaumann derivative, the continuum under a constant shear strain rate $\dot{\gamma}$ rotates endlessly at the spin $\dot{\theta}^W = -\frac{\dot{\gamma}}{2}$. Consequently, the spin of directors is equal to $\dot{\theta}_{dir}^W = \dot{\theta}^W = -\frac{\dot{\gamma}}{2}$. The value of $\sigma_{12}$ oscillates between $\frac{R_0}{\sqrt{2N}}$ at $F_{12} = \gamma = k\pi$ ($k \in \mathbb{N}$) and a peak (or valley) value
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Tension/compression test for two deformation rates $\dot{F}_{11} = 0.1 \ s^{-1}$ (a) and $\dot{F}_{11} = 0.001 \ s^{-1}$ (b). Model with a time-independent and a viscoplastic mechanism</figcaption></figure>

$R_0 / \sqrt{F+G+4H}$ at $\gamma = \frac{2k+1}{2}\pi$ ($k \in \mathbb{N}$). In the case of the Green–Naghdi formulation, the spin of directors is equal to $\dot{\theta}_{dir}^R = -\frac{2\dot{\gamma}}{4+\gamma^2}$. Therefore, the rotation angle converges to $\pi/2$ as $\gamma$ goes to infinity. This is illustrated by Fig. 14.

The responses of the proposed model are shown in Fig. 15. If the plastic spin vanishes i.e. $\beta = 0$, for small elastic stretches, the proposed model without plastic spin and the model based on Jaumann derivative have the same response. This means that the orthotropic axes rotate at the same rate as the continuum $\dot{\theta} = -\dot{\gamma}/2$. In contrast, as shown in Fig. 14, when $\beta \neq 0$MPa$^{-1}$, the spin of the substructure is different from material spin. In fact, as the value of $\beta$ increases, the spin induced by $\overline{W}^p$ balances out the material spin. Further, if $\beta$ is sufficiently high, the rotation of directors saturates rapidly at an angle of $\approx \pi/4$. The rate of directors' spin is illustrated in Fig. 15 for several values of the plastic spin parameter $\beta$: constant spin for $\beta = 0$, rapid saturation for high values and oscillatory response for intermediate values.

## 4.6 Plastic spin: crystal plasticity

Crystal plasticity represents one of the few physical situations for which the plastic spin of crystal directors is precisely known. The spin of directors uniquely results from the slip of $N^s$ systems on specific crystallographic planes and along specific slip directions. The dissipation potential is expressed in terms of the Schmid yield function $f^s$ associated with each slip system $s$

$$f^s = |\tilde{M} : \tilde{N}^s - x^s| - \tau_c^s \quad \text{with} \quad \tilde{N}^s = \underline{I}^s \otimes \tilde{n}^s \quad (96)$$

where $\tau_c^s$ denotes the critical resolved shear stress (CRSS) for the s-th slip system, $\underline{I}^s$ and $\tilde{n}^s$ are, respectively, the slip direction and the normal to the slip plane. The resolved Mandel shear stress $\tau^s = \tilde{M} : \tilde{N}^s$ is the driving force for activation of the s-th slip system. Kinematic hardening has been introduced by [78] in the crystal plasticity framework, in the form of a back stress variable $x^s$ obeying the following evolution rule

$$\dot{\alpha}^s = \dot{\gamma}^s - D|\dot{\gamma}^s|\alpha^s, \quad x^s = C\alpha^s \quad (97)$$

where $C, D$ are the kinematic hardening material parameters. For the sake of demonstration, a power law potential is considered

$$\Omega(\tilde{M}, \tilde{N}^s) = \sum_{s=0}^{N_s} \frac{K}{n+1} \left\langle \frac{|\tilde{M} : \tilde{N}^s - x^s| - \tau_c^s}{K} \right\rangle^{n+1} \quad (98)$$

where $n$ and $K$ are viscosity material parameters and $< \cdot >$ denotes the Macauley brackets. Further, a non-linear hardening law is adopted for the CRSS $\tau_c^s$ given by

$$\tau_c^s = \tau_c^{\mathcal{F}(s)} + \sum_{r=0}^{N_s} H^{sr} (1 - \exp(-b^{\mathcal{F}(r)} \gamma_{cum}^r)) \quad (99)$$

$\mathcal{F}(r)$ identifies the slip system family to which the slip system $r$ belongs (for example basal and prismatic system families in HCP crystals), $N_s$ is the total number of slip systems and $\gamma_{cum}^r$ denotes the accumulated plastic slip. The matrix $H^{sr}$ accounts for interactions between slip systems. The evolution of the plastic slip variables $\gamma^s$ follows as

$$\dot{\gamma}^s = \left\langle \frac{|\tilde{M} : \tilde{N}^s - x^s| - \tau_c^s}{K} \right\rangle^n \text{sign}(\tilde{M} : \tilde{N}^s - x^s), \\ \dot{\gamma}_{cum}^s = |\dot{\gamma}^s| \quad (100)$$
<figure><img src="image_14.png" /><figcaption>Fig. 13 Hypoelastoplastic material response under simple glide loading. Jaumann (left) and Green–Naghdi (right) formulations. Coefficients of Hill’s yield criterion used in this simulation: $F = G = H = 0.5$, $L = M = 1.5$, $N = 4$, initial yield stress $R_0 = 1000$ MPa</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14 Effect of plastic spin on stress values in the case of simple glide loading. a Without plastic spin i.e. $\beta = 0$ MPa$^{-1}$ b with plastic spin $\beta = 0.01$ MPa$^{-1}$. Multiplicative plasticity model with yield stress $R_0 = 1000$ MPa, Hill yield criterion coefficients $F = G = H = 0.5$, $L = M = 1.5$, $N = 4$</figcaption></figure>

According to the normality rule (27), the plastic deformation rate reads

$$ \overline{L}^p = \sum_{s=0}^{N_s} \dot{\gamma}^s \underline{N}^s \qquad (101) $$

and the plastic spin writes

$$ \overline{W}^p = \sum_{s=0}^{N_s} \dot{\gamma}^s \text{skw} (\underline{l}^s \otimes \underline{n}^s) \qquad (102) $$

As an illustration, the simple glide kinematics (84) is imposed to a face centered cubic (FCC) single crystal. The single crystal model response is compared to two fake crystal models in order to highlight the importance of properly characterizing the plastic spin. In the fake models no plastic spin is introduced and Jaumann or Green–Naghdi hypoelastic laws are used instead of the actual plastic spin (102). The responses

<figure><img src="image_16.png" /><figcaption>Fig. 15 Directors’ spin for various values of parameter $\beta$ during a simple glide test for a Hill perfectly plastic material and multiplicative decomposition. When $\beta = 0$MPa$^{-1}$, the rotation rate of directors is constant and coincides with the material spin</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16 The effect of plastic spin in the case of a simple glide loading of a FCC single crystal for three different finite strain formulations: Jaumann (a) and Green–Naghdi (b) rates (with no plastic spin) and the multiplicative decomposition (c). (d) comparison of $\sigma_{12}$ for all three formulations</figcaption></figure>

of the three models are shown in Fig. 16. The material is perfectly plastic, with an initial CRSS value $\tau_c^s = 100$ MPa. The simple glide loading direction 1 initially coincides with the crystal direction [100], whereas the direction 2 initially coincides with the crystal direction [0$\bar{1}$1]. The octahedral slip system family of the FCC crystal is employed and contains 12 slip systems. The reference solution based on the multiplicative decomposition and plastic spin predicts a saturation of stress levels after $F_{12} = 2$ shear. In contrast, the fake crystal responses exhibit oscillatory stress evolutions. It is noteworthy that the stress component $\sigma_{33}$ vanishes for hypoelastic models. On the other hand, the approach based on the multiplicative decomposition includes a plastic spin derived from the dissipation potential. This model of plastic spin dictates the spin of directors independently of material rotation (see Fig. 17). In addition, the $\sigma_{33}$ does not vanish according to the present model contrary to the hypoelastic formulations.

Furthermore, the rotation rate tends to zero for the crystal orientation considered in the present example, which is not the case for Jaumann rate. The comparison with experimental results confirms that the rotation of anisotropic axes

of a single crystal does not follow the material rotation [79, 80]. Clearly, the fake crystal plasticity models are physically inadequate. Note that oscillatory responses can be observed for special crystal orientations within the framework of multiplicative crystal plasticity, see [81].

## 4.7 Crystal plasticity with system interactions

The interaction matrix in Eq. (99) has the general form

$$ H = \begin{pmatrix} Q_1 H_1 & Q_1 h_2^1 & \dots & Q_1 h_n^1 \\ Q_2 h_2^1 & Q_2 H_2 & \dots & Q_2 h_n^2 \\ \vdots & \vdots & \ddots & \vdots \\ Q_n h_n^1 & Q_n h_n^2 & \dots & Q_n H_n \end{pmatrix} \quad (103) $$

where $Q_i$ is the isotropic hardening parameter associated with each slip system family $i$, $H_i$ matrices denote the self-
<figure><img src="image_18.png" /><figcaption>Fig. 17 The spin of material directors in the case of a simple glide loading for three different finite strain formulations of crystal plasticity: Jaumann and Green–Naghdi rates and the multiplicative decomposition</figcaption></figure>

hardening coefficients, and $h_m^n$ denotes the latent hardening parameters (hardening of slip systems belonging to the family $n$ caused by slip systems belonging to the family $m$). The table 2 shows typical material constants for a single crystal with octahedral and cubic slip system families. Such combination of octahedral and cube slip system families are encountered in single crystal Nickel–based superalloys [78].

Figure 18 shows the influence of latent hardening parameter $h_2^1$ on the activation of slip systems. In fact, if $h_2^1 = 0$, the cubic slip systems are activated as soon as $\sigma_{33} \approx \frac{3\sqrt{2}}{2}\tau_c^{(2)}$ (according to Schmid’s law, neglecting the overstress due to viscoplasticity), and the octahedral slip systems are activated when $\sigma_{33} \approx \frac{3\sqrt{6}}{2}\tau_c^{(1)}$. If $h_2^1 > 0$, the octahedral slip systems will be activated later since the corresponding critical resolved shear increases due to accumulated cubic plastic slip according to Eq. (99). This is the manifestation of latent hardening between slip system families.

## 5 Structural applications

The proposed generic formulation and implementation are now illustrated in the case of structural components subjected to various loading conditions. The model predictions and computational efficiency are compared to those obtained by

standard approaches involving hypo-elastoviscoplastic models.

### 5.1 Application 1: Deep drawing for anisotropic materials

The present approach is applied first to the three dimensional problem of cup deep drawing. This problem is common in literature and solved using various finite strain formulations, see [23, 44]. The geometry of the test is described in Fig. 19.

A time-independent elastoplastic model is considered with a Hill yield function as in Eq. (95). A nonlinear isotropic hardening rule is adopted in the form

$$R(p) = R_0 + Q(1 - \exp(-bp)) \quad (104)$$

Here the value of the initial yield stress is $R_0 = 253$ MPa, the isotropic hardening parameters are chosen as $b = 14$ and $Q = 215$ MPa. The parameters of Hill’s yield function are set to $F = G = H = N = 0.5$ and $L = M = 1.5$.

Herein, only a quarter of the cup is analyzed due to the orthotropic material symmetry. The mesh of the sheet contains 900 C3D8 hexahedral solid elements (linear interpolation with 8 Gauss points per element), with 3 elements in the thickness direction. The plate is initially circular. The tools are modeled as rigid surfaces and are completely fixed except the punch which is pushed in z-direction to a total displacement of $u_z = -40$ mm.

Figure 20 shows the distribution of the accumulated plastic strain and the deformed configuration. If the plastic spin vanishes, the maximal values of the accumulated plastic strain are located along the directions of the material symmetry x- and y- axes. When plastic spin is taken into account by means of parameter $\beta$ from Eq. (45), the distribution of accumulated plastic strain tends to be more isotropic for $\beta = 0.01$ MPa$^{-1}$.

The earing profile after forming is depicted in Fig. 21 for two values of the parameter $\beta$. The edge of the sheet has a wavy shape. As the value of $\beta$ increases, the edge shows less pronounced wave-shape. Similar results were reported in [44] for an isotropic elastic and a transverse isotropic Hill-type plastic model. The hypoelastic version of this model yields similar results to the present approach, including the wavy shape of the cup and the distribution of plastic strains. The present formulation, however, is advantageous since a supplementary plastic spin describing the evolution of anisotropy axes can be incorporated in the model.

<table><caption>Table 2 Material parameters for octahedral (1) and cubic (2) slip systems with isotropic hardening</caption><thead><tr><th>K<sub>1</sub><br>(MPa<sup>1/n<sub>1</sub></sup>)</th><th>n<sub>1</sub><br>(–)</th><th>τ<sub>c</sub><sup>(1)</sup><br>(MPa)</th><th>Q<sub>1</sub><br>(MPa)</th><th>b<sub>1</sub><br>(–)</th><th>H<sub>1</sub><br>(–)</th><th>K<sub>2</sub><br>(MPa<sup>1/n<sub>2</sub></sup>)</th><th>n<sub>2</sub><br>(–)</th><th>τ<sub>c</sub><sup>(2)</sup><br>(MPa)</th><th>Q<sub>2</sub><br>(MPa)</th><th>b<sub>2</sub><br>(–)</th><th>H<sub>2</sub><br>(–)</th></tr></thead><tbody><tr><td>0.1</td><td>20</td><td>100</td><td>50</td><td>50</td><td>1</td><td>0.1</td><td>20</td><td>100</td><td>100</td><td>100</td><td>1</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>**Fig. 18** A single crystal including two slip system families (octahedral an cubic) in tension along < 111 >. $\gamma_{cum}^{cub}$ (resp. $\gamma_{cum}^{oct}$) denotes the accumulated cubic (resp. octahedral) plastic slip</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19** Schematic view of the cup drawing process. No friction is considered between sheet and tools. The die, punch and the holder are assumed to be linear elastic with a Young modulus of $E = 10^6$ MPa and a Poisson ratio $\nu = 0.3$. Lengths in mm</figcaption></figure>

## 5.2 Application 2: Turbine blade with single crystal plasticity

In the following, the behavior of a nickel-based superalloy single crystal turbine blade subjected to creep at high temperature is studied. During their operation, turbine blades are subjected to centrifugal forces induced by the rotation (~ 20000 RPM) of the turbine disc in addition to gas pressure. During one flight, the turbine blades are subjected to high and non-uniform temperatures (maximum temperature ~ 1200 °C), which will induce thermal strains. The thermomechanical loading is maintained for a longer time, compared to the nominal in-service conditions, at the maximum temperature(~ 640 °C at the root and ~ 1200 °C at the tip). The mesh of the turbine blade geometry contains 1366 linear C3D8 hexahedral solid elements. We consider two initial orientations of the crystal. The first one is such that the crystal directions triplet ([100]–[010]–[001]) coincide with the orthogonal basis vectors triplet ($x_1,x_2,x_3$). For the second orientation, the crystal is tilt by 15° in the $x_1-x_3$ plane. The

objective is to assess the impact of crystal misorientation on the thermomechanical response of the blade.

The material model has 18 slip systems potentially active, 12 octahedral {110} <111> and 6 cubic slip systems {110} <100>. The constitutive equations of the model including kinematic hardening law were given in Sect. 4.6. Typical values of the material parameters used in the simulation can be found in [82]. Cubic elasticity moduli, the coefficient of isotropic thermal expansion, critical resolved shear stresses, viscoplastic flow parameters and nonlinear kinematic hardening parameters are identified as functions of temperature.

Figure 22 shows that both hypoelastic formulations yield similar creep results. The present model prediction slightly differs and the difference increases with time. The difference between various formulations becomes apparent when an initial rotation of 15° around y-axis is considered. Figure 23 shows the relative rotation of material directors with respect to material rotation for the present model. For hypoelastic formulations, the spin of directors coincides with material rotation which is not the case for the model based on the multiplicative decomposition. The multiplicative crystal plasticity model assumes a relative rotation of crystal directors with respect to material lines, induced by the plastic spin (102). An initial misorientation of 15° leads to slightly larger relative rotation of crystal directors, as shown in Fig. 23. This result explains the fact that the difference in creep predicted by various formulations is much more pronounced for an initial misorientation of 15° compared to 0°. It must be underlined that the presented crystal plasticity formulations based on corotational frames are physically unsound but represent standard extensions of small strain crystal plasticity.
<figure><img src="image_21.png" /><figcaption>**Fig. 20** Accumulated plastic strain induced by cup drawing in an anisotropic elastoplastic plate: **a** Without plastic spin $\beta = 0 \text{ MPa}^{-1}$, **b** with plastic spin $\beta = 0.01 \text{ MPa}^{-1}$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 21** Distance of rim to the sheet center versus angle to x-axis. The curve shows the earing pattern for two values of the parameter $\beta$</figcaption></figure>

## 5.3 Computational efficiency of the approach

The present formulation is compared to two hypoelastic formulations, Jaumann (J) and Green–Naghdi (GN) rates, in terms of computational efficiency. We consider four examples of structural applications and different material models. In addition to the applications discussed previously, we carry out a tensile test on a notched specimen with 323703 nodes corresponding to 73920 C3D20R quadratic elements with reduced integration. A von Mises plasticity model is adopted with initial yield stress $R_0 = 300 \text{ MPa}$. The parameter values for nonlinear isotropic and kinematic hardening are taken as $Q = 400 \text{ MPa}$, $b = 2.5$, $C = 5700 \text{ MPa}$ and $D = 17$, as in [33]. The results predicted by different formulations (hypoelastic and multiplicative models) turn out to be almost identical in this isotropic case. Figure 24 depicts the obtained accumulated plastic strain field. The objective of this example

<figure><img src="image_23.png" /><figcaption>**Fig. 22** **a** The accumulated plastic strain in a turbine blade subjected to creep at high temperature. **b** Displacement of the node Q as a function of time: comparison of different finite deformation formulations</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 23** Relative rotation of crystal directors w.r.t. material lines for the multiplicative-based model in two cases a ideal < 1001 > orientation of the blade, b deviation of 15° around y-axis</figcaption></figure>

is to illustrate the performance of the approach for a problem with almost one million degrees of freedom.

Note that the integration of constitutive equations is carried out using a fully implicit integration method ($\theta = 1$). The time increment is set to the same value for all the formulations (if it does not converge, the time increment is divided by 2). Comparison results for various formulations are presented in Fig. 25. Computations were carried out on processors of type Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz. It is found that the formulation based on the multiplicative decomposition leads to a slightly higher computational cost compared to hypoelastic formulations (+15–20%). This is due to the fact that more operations are implied by the systematic use of non-symmetric tensors in the algorithm. Besides, the total number of iterations for convergence of global equilibrium, is slightly higher for the model based on the multiplicative decomposition. In spite of that, the computational cost is still reasonable compared to formulations common in commercial FEM software. The slight increase in computation time is counter-balanced by additional possibilities in the modeling of anisotropic inelasticity.

# 6 Conclusion

The present work demonstrates that elasto-viscoplastic models based on the multiplicative decomposition of the deformation gradient are now mature for a systematic use in commercial finite element codes for structural computations.

<figure><img src="image_25.png" /><figcaption>**Fig. 24** Accumulated plastic strain in a notched specimen subjected to tensile loading using the multiplicative decomposition based model</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 25** Comparison of the total number of increments, total number of iterations and CPU time for different formulations. Applications: a deep drawing with 900 C3D8 elements (1324 nodes), **b** deep drawing with</figcaption></figure>

19500 C3D8 elements (26564 nodes), **c** turbine blade under creep with 1366 C3D8 elements (2498 nodes), **d** Notched specimen with 73920 C3D20R elements (323703 nodes). Simulations are run on 24 CPUs

The proposed generic constitutive framework overcomes the shortcomings of standard formulations used in most available FE codes and based on hypoelastic laws and limited description of anisotropic behavior. It has been shown that the computing efficiency is comparable to the standard approach, although slightly less advantageous. The use of a thermodynamically consistent formulation of constitutive equations ensures increased reliability of model predictions. It also contains new model formulations including multiple inelastic mechanisms and internal hardening variable that remain to be applied to available experimental data.

The proposed extension of constitutive equations at finite strain relies upon the multiplicative decomposition of the deformation gradient and a hyperelastic relation between stress and elastic strain measures. The other characteristic feature is the expression of plasticity laws in terms of the Mandel stress tensor, which is generally non-symmetric. The approach accommodates contributions of many deformation mechanisms combining thermo-plasticity, viscosity and possibly damage, without resorting to arbitrary further multiplicative decompositions of elastic or plastic parts of the deformation gradient. Instead, the rates of multimechanism contributions are added in the viscoplastic flow rule. Anisotropy is incorporated via the consideration of directors and associated structural tensors and plastic spin concept

which are often absent in existing standard formulations. It was illustrated in the case of Hill's criterion and crystal plasticity.

Two numerical integration methods of the nonlinear evolution equations are proposed. First, the elastic or the plastic part of the deformation gradient is integrated using an exponential map. In general, this method involves calculating the exponential of a non-symmetric second-order tensor and its derivative, which is not an easy task. An alternative integration method is considered, which consists of correcting the elastic part of the deformation gradient to enforce the plastic incompressibility. The present framework is illustrated through various models, including isotropic and anisotropic (visco)-plasticity with isotropic and kinematic hardening. Most formulations exist in the literature but a few are original: nonlinear non-symmetric kinematic hardening rule including dynamic and static recovery terms. In the case of linear kinematic hardening, models based on the Jaumann rate and the multiplicative decomposition exhibit stress oscillations extensively discussed in the literature. However, regarding nonlinear kinematic hardening including dynamic recovery term, the various models provide similar predictions for a proper choice of material parameters.

The capabilities of the present formulation are illustrated through elementary industrial applications. The effect of
plastic spin on the evolution of anisotropy described by a triad of directors was evaluated in the case of cup drawing. This effect cannot be neglected, particularly for materials showing a high degree of anisotropy under large plastic deformation. The comparison between hypo-elastoplastic models applied to crystal plasticity shows that the result heavily depends on the objective stress rate. The crystal plasticity model based on the multiplicative decomposition remains the reference model in that case and was successfully validated by experiments, especially the rotation of crystal directors with respect to loading axes. The model formulation should also be applied to benchmarks available in finite elastoplasticity as proposed in [83]. Consideration of enhanced element formulations such mixed or hybrid elements [84] is recommended to improve the treatment of plastic incompressibility.

Meanwhile, further research is necessary to develop yield criteria in terms of the generally non-symmetric Mandel stress tensor [44, 46]. The identification of the constitutive equations for the plastic spin requires more experimental investigations at large deformations. Furthermore, non-symmetric internal variables (e.g. kinematic hardening) are to be considered. Examples showing applications to compressible elastoviscoplastic materials could also be provided, for instance based on the family of yield functions proposed and implemented in [85]. Non-associated flow rules for applications to granular media should also be envisaged.

It is hoped that the present results can increase the interest of the engineering computation community towards the systematic use of physically consistent nonlinear constitutive equations and improved modeling of anisotropic plasticity at finite deformation.

## Appendix A Hypoelastic-based formulations

The corotational strain rate $\boldsymbol{\varepsilon}$ and stress $\boldsymbol{s}$ are given by

$$ \dot{\boldsymbol{\varepsilon}} = \boldsymbol{\varepsilon}^T \boldsymbol{D} \boldsymbol{\varepsilon}, \quad \boldsymbol{s} = \boldsymbol{\varepsilon}^T \boldsymbol{\sigma} \boldsymbol{\varepsilon} \quad (A1) $$

where $\boldsymbol{Q}$ is a rotation tensor. For Jaumann-rate based formulation, $\boldsymbol{\tilde{Q}}$ is obtained from the integration of the differential equation $\boldsymbol{W} = \dot{\boldsymbol{Q}} \boldsymbol{\tilde{Q}}^T$ where $\boldsymbol{W} = \text{skw}(\boldsymbol{L})$ is the skew-symmetric part of the velocity gradient. For Green-Naghdi rate, $\boldsymbol{Q} = \boldsymbol{R}$, where $\boldsymbol{R}$ is an orthogonal tensor obtained from the polar decomposition $\boldsymbol{F} = \boldsymbol{R}\boldsymbol{U}$ ($\boldsymbol{U}$ is the right stretch tensor). Accordingly, the spin tensors corresponding to Jaumann and Green-Naghdi rates are defined as

$$ \boldsymbol{\Omega}^J = \boldsymbol{W} \quad \text{and} \quad \boldsymbol{\Omega}^{GN} = \dot{\boldsymbol{R}} \boldsymbol{R}^T \quad (A2) $$

The corotational strain rate is additively split into elastic and plastic parts:

$$ \dot{\boldsymbol{\varepsilon}} = \dot{\boldsymbol{\varepsilon}}^e + \dot{\boldsymbol{\varepsilon}}^p \quad (A3) $$

The elasticity law is defined as

$$ \boldsymbol{s} = \mathcal{C}_{\tilde{\sim}} : \boldsymbol{\varepsilon}^e \quad (A4) $$

where $\mathcal{C}_{\tilde{\sim}}$ is a fourth-order elasticity moduli. The elasticity law can be written in rate form, as an hypoelastic constitutive law involving an objective derivative:

$$ \frac{D\boldsymbol{\sigma}}{Dt} = \dot{\boldsymbol{\varepsilon}} + \boldsymbol{\sigma}\boldsymbol{\Omega} - \boldsymbol{\Omega}\boldsymbol{\sigma}, \quad \text{with} \quad \boldsymbol{\Omega} = \boldsymbol{\Omega}^J, \boldsymbol{\Omega}^{GN} \quad (A5) $$

corresponding respectively to the Jaumann and Green-Naghdi rates. The yield condition is given by

$$ f(\boldsymbol{s}, \boldsymbol{A}_I, T) = 0 \quad (A6) $$

The plastic part of the corotational strain rate is determined by the flow rule

$$ \dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial f(\boldsymbol{s}, \boldsymbol{A}_I, T)}{\partial \boldsymbol{s}} \quad (A7) $$

## Appendix B Finite element implementation

The nodal displacements are written in vector form as

$$ \{\boldsymbol{u}^e\} = \{u_1^1 \ u_2^1 \ u_3^1 \ \dots \ u_1^p \ u_2^p \ u_3^p\}^T \quad (B8) $$

$u_j^i$ denotes the j-th ($1 \le j \le 3$) component of the nodal displacement at node i. $p$ is the total number of nodes. The shape functions can be written as

$$ [N] = \begin{bmatrix} N^1 & 0 & 0 & \dots & N^p & 0 & 0 \\ 0 & N^1 & 0 & \dots & 0 & N^p & 0 \\ 0 & 0 & N^1 & \dots & 0 & 0 & N^p \end{bmatrix} \quad (B9) $$

and

$$ \{\boldsymbol{u}\} = [N]\{\boldsymbol{u}^e\} \quad (B10) $$

It follows that

$$ \begin{aligned} \{\boldsymbol{F}\} &= \{F_{11} \ F_{22} \ F_{33} \ F_{12} \ F_{23} \ F_{31} \ F_{21} \ F_{32} \ F_{13}\}^T \\ &= [\boldsymbol{B}][\boldsymbol{u}^e] \end{aligned} \quad (B11) $$
where

$$[\mathbf{B}] = \begin{bmatrix} \frac{\partial N_1}{\partial X_1} & 0 & 0 & \dots & \frac{\partial N^p}{\partial X_1} & 0 & 0 \\ 0 & \frac{\partial N_1}{\partial X_2} & 0 & \dots & 0 & \frac{\partial N^p}{\partial X_2} & 0 \\ 0 & 0 & \frac{\partial N_1}{\partial X_3} & \dots & 0 & 0 & \frac{\partial N^p}{\partial X_3} \\ \frac{\partial N_1}{\partial X_2} & 0 & 0 & \dots & \frac{\partial N^p}{\partial X_2} & 0 & 0 \\ 0 & \frac{\partial N_1}{\partial X_3} & 0 & \dots & 0 & \frac{\partial N^p}{\partial X_3} & 0 \\ 0 & 0 & \frac{\partial N_1}{\partial X_1} & \dots & 0 & \frac{\partial N^p}{\partial X_1} & 0 \\ 0 & \frac{\partial N_1}{\partial X_1} & 0 & \dots & 0 & \frac{\partial N^p}{\partial X_1} & 0 \\ \frac{\partial N_1}{\partial X_3} & 0 & 0 & \dots & \frac{\partial N^p}{\partial X_3} & 0 & 0 \end{bmatrix} \quad (\text{B12})$$

The linear set of equations to be solved for $\Delta \underline{\mathbf{u}}$ iteratively is given by

$$[\mathbf{K}]^g \Delta \underline{\mathbf{u}} = \int_{\Omega_0} [N^T]\{\mathbf{b}\} dV_0 + \int_{\partial \Omega_0^T} [N]^T\{\mathbf{t}\} dS_0 \\ - \int_{\Omega_0} [\mathbf{B}]^T\{\mathbf{P}\} dV_0 \quad (\text{B13})$$

where $[\mathbf{K}]^g$ is the global tangent matrix in Eq. (63) and Eq. (67) and $\{\mathbf{P}\}$ is the Boussinesq stress tensor. $\{\mathbf{b}\}$ and $\{\mathbf{t}\}$ denote the reference body force and surface traction fields, respectively.

## Appendix C Tangent matrices

The Jacobian matrix is required to integrate the constitutive equations at the Gauss point level. The block form of the Jacobian matrix writes

$$\begin{pmatrix} \frac{\partial \mathcal{R}_{el}}{\partial \Delta \mathbf{F}^e} & \frac{\partial \mathcal{R}_{el}}{\partial \Delta p_i} & \frac{\partial \mathcal{R}_{el}}{\partial \Delta \boldsymbol{\alpha}_i} \\ \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \mathbf{F}^e} & \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta p_i} & \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \boldsymbol{\alpha}_i} \\ \frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta \mathbf{F}^e} & \frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta p_i} & \frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta \boldsymbol{\alpha}_i} \end{pmatrix} \quad (\text{C14})$$

where the residuals are taken as

$$\mathcal{R}_{el} = \Delta \mathcal{F}^e - \Delta \mathcal{F} \mathcal{F}^{-1} \mathcal{F}^e + \mathcal{F}^e \sum_{i=0}^{N} \Delta p_i \mathcal{N}_i \quad (\text{C15})$$

$$\mathcal{R}_{p_i} = f_i \quad \text{or} \quad \mathcal{R}_p = \Delta p_i - \Delta t \frac{\partial \Omega}{\partial f_i} \quad (\text{C16})$$

$$\mathcal{R}_{\boldsymbol{\alpha}_i} = \Delta \boldsymbol{\alpha}_i - \Delta p_i \mathbf{m}_i + \Delta \mathbf{q}_i \quad (\text{C17})$$

The terms related to the global part of the Jacobian (the first row and the first column) are given by

$$\frac{\partial \mathcal{R}_{el}}{\partial \Delta \mathbf{F}^e} = \mathbf{1} - (\Delta \mathcal{F} \mathcal{F}^{-1}) \otimes \mathbf{1} + \sum_{i=0}^{N} \Delta p_i (\mathbf{z} \otimes \mathcal{N}_i^T)$$

$$+ \theta \mathcal{F}^e \left[ \sum_{i=0}^{N} \Delta p_i \frac{\partial \mathcal{N}_i}{\partial \mathcal{M}} \right] : \frac{\partial \mathcal{M}}{\partial \mathcal{F}^e} \quad (\text{C18})$$

$$\frac{\partial \mathcal{R}_{el}}{\partial \Delta p_i} = \mathcal{F}^e \mathcal{N}_i + \Delta p_i \frac{\partial \mathcal{N}_i}{\partial \mathcal{A}_i} \frac{\partial \mathcal{A}_i}{\partial p_i} \quad (\text{C19})$$

$$\frac{\partial \mathcal{R}_{el}}{\partial \Delta \boldsymbol{\alpha}_i} = \theta \mathcal{F}^e \sum_{i=0}^{N} \Delta p_i \left( \frac{\partial \mathcal{N}_i}{\partial \mathcal{A}_i} \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_i} \right) \quad (\text{C20})$$

$$\frac{\partial \mathcal{R}_p}{\partial \Delta \mathcal{F}^e} = \theta \frac{\partial f}{\partial \mathcal{M}} : \frac{\partial \mathcal{M}}{\partial \mathcal{F}^e} \quad \text{or}$$

$$\frac{\partial \mathcal{R}_p}{\partial \Delta \mathcal{F}^e} = -\theta \Delta t \frac{\partial \dot{p}_i}{\partial f_i} \left( \frac{\partial f}{\partial \mathcal{M}} : \frac{\partial \mathcal{M}}{\partial \mathcal{F}^e} \right) \quad (\text{C21})$$

$$\frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta \mathcal{F}^e} = -\theta \Delta p_i \frac{\partial \mathbf{m}_i}{\partial \mathcal{M}} : \frac{\partial \mathcal{M}}{\partial \mathcal{F}^e} \quad (\text{C22})$$

where

$$\frac{\partial \mathcal{M}}{\partial \mathcal{F}^e} = \mathbf{1} \otimes (\mathbf{\Pi}^e \mathcal{F}^e T) + \mathcal{F}^e T \otimes \mathbf{\Pi}^e + \mathcal{F}^e T \mathcal{F}^e \frac{\partial \mathbf{\Pi}^e}{\partial \mathcal{F}^e} : \frac{\partial \mathcal{F}^e}{\partial \mathcal{F}^e} \quad (\text{C23})$$

$$\frac{\partial \mathcal{F}^e}{\partial \mathcal{F}^e} = \frac{1}{2} \left( \mathbf{1} \otimes \mathcal{F}^e T + \mathcal{F}^e T \otimes \mathbf{1} \right) \quad (\text{C24})$$

Recall that $\theta \in [0,1]$ is a parameter of the integration method. Next, the block of the Jacobian related to the internal variable evolution equations writes

$$\frac{\partial \mathcal{R}_{p_i}}{\partial \Delta p_i} = \theta \frac{\partial f}{\partial p} \quad \text{or} \quad \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta p_i} = 1 - \theta \Delta t \frac{\partial \dot{p}_i}{\partial f_i} \frac{\partial f_i}{\partial p_i} \quad (\text{C25})$$

$$\frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \boldsymbol{\alpha}_i} = \frac{\partial f}{\partial \boldsymbol{\alpha}_i} \frac{\partial \boldsymbol{\alpha}_i}{\partial \Delta \boldsymbol{\alpha}_i} \quad \text{or} \quad \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \boldsymbol{\alpha}_i} = -\theta \Delta t \frac{\partial \dot{p}_i}{\partial f_i} \frac{\partial f_i}{\partial \mathcal{A}_i} \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_i} \quad (\text{C26})$$

$$\frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta p_i} = -\mathbf{m}_i - \theta \Delta p_i \frac{\partial \mathbf{m}_i}{\partial p_i} \quad (\text{C27})$$

$$\frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta \boldsymbol{\alpha}_i} = \mathbf{1} - \theta \left( \Delta p_i \frac{\partial \mathbf{m}_i}{\partial \mathcal{A}_i} + \Delta t \frac{\partial \mathbf{q}_i}{\partial \mathcal{A}_i} \right) \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_i} - \theta \Delta p_i \frac{\partial \mathbf{m}_i}{\partial \boldsymbol{\alpha}_i} \quad (\text{C28})$$

The interaction terms are given by

$$\frac{\partial \mathcal{R}_{p_i}}{\partial \Delta p_j} = \theta \frac{\partial f_i}{\partial p_j} \quad \text{or} \quad \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta p_j} = -\theta \Delta t \frac{\partial \dot{p}_i}{\partial f_i} \frac{\partial f_i}{\partial p_j} \quad (\text{C29})$$

$$\frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \boldsymbol{\alpha}_i} = \theta \frac{\partial f}{\partial \mathcal{A}_i} \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_j} \quad \text{or} \quad \frac{\partial \mathcal{R}_{p_i}}{\partial \Delta \boldsymbol{\alpha}_i} = -\theta \Delta t \frac{\partial \dot{p}_i}{\partial f_i} \frac{\partial f_i}{\partial \mathcal{A}_i} \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_j} \quad (\text{C30})$$

$$\frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta p_i} = \mathbf{0} \quad (\text{C31})$$

$$\frac{\partial \mathcal{R}_{\boldsymbol{\alpha}_i}}{\partial \Delta \boldsymbol{\alpha}_j} = \mathbf{1} - \theta \left( \Delta p_i \frac{\partial \mathbf{m}_i}{\partial \mathcal{A}_i} + \Delta t \frac{\partial \mathbf{q}_i}{\partial \mathcal{A}_i} \right) \frac{\partial \mathcal{A}_i}{\partial \boldsymbol{\alpha}_j} \quad (\text{C32})$$
The matrix $\frac{\partial A_i}{\partial \alpha_j}$ accounts for interactions between the hardening variables.

The second method to integrate constitutive equations relies upon the following definition of the residual in Eq. (C15)

$$ \mathcal{R}_{el} = \mathcal{F}^e \mathcal{F}^p - \mathcal{F} \qquad (C33) $$

$\mathcal{F}^p$ is the solution of the following differential equation

$$ \dot{\mathcal{F}}^p = \overline{\mathcal{L}}^p \mathcal{F}^p \qquad (C34) $$

So $\mathcal{F}^p$, at increment $n + 1$, can be estimated by using the exponential mapping [75] as

$$ \mathcal{F}_{n+1}^p = \exp(\Delta \overline{\mathcal{L}}^p) \mathcal{F}_n^p \qquad (C35) $$

Accordingly, the plastic incompressibility, i.e. trace($\overline{\mathcal{L}}^p$) = 0 or det($\mathcal{F}^p$) = 1, is satisfied since

$$ \det(\exp \mathcal{A}) = \exp(\operatorname{tr} \mathcal{A}) \qquad (C36) $$

The first row of the Jacobian matrix in Eq. (C14) is rewritten

$$ \frac{\partial \mathcal{R}_{el}}{\partial \Delta \mathcal{F}^e} = \mathbf{1} \otimes \mathcal{F}^p \mathbf{T} + \theta \mathcal{F}^e \left[ \frac{\partial \exp(\Delta \overline{\mathcal{L}}^p)}{\partial \Delta \overline{\mathcal{L}}^p} : \left( \sum_{i=0}^{N} \Delta p_i \frac{\partial \mathcal{N}_i}{\partial \overline{\mathcal{M}}} \right) : \frac{\partial \overline{\mathcal{M}}}{\partial \mathcal{F}^e} \right] \mathcal{F}_n^p \qquad (C37) $$

$$ \frac{\partial \mathcal{R}_{el}}{\partial \Delta p_i} = \mathcal{F}^e \left( \frac{\partial \exp(\Delta \overline{\mathcal{L}}^p)}{\partial \Delta \overline{\mathcal{L}}^p} : \mathcal{N}_i \right) \mathcal{F}_n^p \qquad (C38) $$

$$ \frac{\partial \mathcal{R}_{el}}{\partial \Delta \alpha_i} = \mathcal{F}^e \left[ \frac{\partial \exp(\Delta \overline{\mathcal{L}}^p)}{\partial \Delta \overline{\mathcal{L}}^p} : \left( \sum_{i=0}^{N} \Delta p_i \frac{\partial \mathcal{N}_i}{\partial \Delta \alpha_i} \right) \right] \mathcal{F}_n^p \qquad (C39) $$

**Funding** Open access funding provided by Mines Paris - PSL.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Green AE, Naghdi PM (1965) A general theory of an elastic-plastic continuum. Arch Ration Mech Anal 18(4):251–281. https://doi.org/10.1007/BF00251666

2. Naghdi PM (1990) A critical review of the state of finite plasticity. Zeitschrift für angewandte Mathematik und Physik ZAMP 41(3):315–394. https://doi.org/10.1007/BF00959986

3. Truesdell C (1955) Hypo-elasticity. J Ration Mech Anal 4:83–1020

4. Hibbitt HD, Marcal PV, Rice JR (1970) A finite element formulation for problems of large strain and large displacement. Int J Solids Struct 6(8):1069–1086. https://doi.org/10.1016/0020-7683(70)90048-X

5. Kojić M, Bathe K-J (1987) Studies of finite element procedures–Stress solution of a closed elastic strain path with stretching and shearing using the updated Lagrangian Jaumann formulation. Comput Struct 26(1):175–179. https://doi.org/10.1016/0045-7949(87)90247-1

6. Lin RC, Schomburg U, Kletschkowski T (2003) Analytical stress solutions of a closed deformation path with stretching and shearing using the hypoelastic formulations. Eur J Mech A Solids 22(3):443–461. https://doi.org/10.1016/S0997-7538(03)00031-7

7. Brepols T, Vladimirov IN, Reese S (2014) Numerical comparison of isotropic hypo- and hyperelastic-based plasticity models with application to industrial forming processes. Int J Plast 63:18–48. https://doi.org/10.1016/j.ijplas.2014.06.003

8. Xiao H, Bruhns OT, Meyers A (1998) On objective corotational rates and their defining spin tensors. Int J Solids Struct 35(30):4001–4014. https://doi.org/10.1016/S0020-7683(97)00267-9

9. Xiao H, Bruhns OT, Meyers A (1999) Existence and uniqueness of the integrable-exactly hypoelastic equation $\dot{\mathcal{L}} = \lambda(\operatorname{tr} \mathcal{D})\mathbf{1} + 2\mu \mathcal{D}$ and its significance to finite inelasticity. Acta Mech 138(1):31–50. https://doi.org/10.1007/BF01179540

10. Jiao Y, Fish J (2017) Is an additive decomposition of a rate of deformation and objective stress rates passé? Comput Methods Appl Mech Eng 327:196–225. https://doi.org/10.1016/j.cma.2017.07.021

11. Dienes JK (1979) On the analysis of rotation and stress rate in deforming bodies. Acta Mech 32(4):217–232. https://doi.org/10.1007/BF01379008

12. Nemat-Nasser S (1979) Decomposition of strain measures and their rates in finite deformation elastoplasticity. Int J Solids Struct 15(2):155–166. https://doi.org/10.1016/0020-7683(79)90019-2

13. Shutov AV, Ihlemann J (2014) Analysis of some basic approaches to finite strain elasto-plasticity in view of reference change. Int J Plast 63:183–197. https://doi.org/10.1016/j.ijplas.2014.07.004

14. Miehe C, Apel N, Lambrecht M (2002) Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials. Comput Methods Appl Mech Eng 191(47):5383–5425. https://doi.org/10.1016/S0045-7825(02)00438-3

15. Friedlein J, Mergheim J, Steinmann P (2022) Observations on additive plasticity in the logarithmic strain space at excessive strains. Int J Solids Struct 239–240:111416. https://doi.org/10.1016/j.ijsolstr.2021.111416

16. Kröner E (1959) Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen. Arch Ration Mech Anal 4(1):273–334. https://doi.org/10.1007/BF00281393

17. Lee EH, Liu DT (1967) Finite strain elastic-plastic theory with application to plane-wave analysis. J Appl Phys 38(1):19–27. https://doi.org/10.1063/1.1708953

18. Mandel J (ed) (1972) Plasticité Classique et Viscoplasticité. Springer, Wien
19. Mandel J (1973) Équations constitutives et directeurs dans les milieux plastiques et viscoplastiques. Int J Solids Struct 9(6):725-740. https://doi.org/10.1016/0020-7683(73)90120-0

20. Lion A (2000) Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. Int J Plast 16(5):469-494. https://doi.org/10.1016/S0749-6419(99)00038-8

21. Wallin M, Ristinmaa M, Ottosen NS (2003) Kinematic hardening in large strain plasticity. Eur J Mech A Solids 22(3):341-356. https://doi.org/10.1016/S0997-7538(03)00026-3

22. Menzel A, Ekh M, Runesson K, Steinmann P (2005) A framework for multiplicative elastoplasticity with kinematic hardening coupled to anisotropic damage. Int J Plast 21(3):397-434. https://doi.org/10.1016/j.ijplas.2003.12.006

23. Vladimirov IN, Pietryga MP, Reese S (2010) Anisotropic finite elastoplasticity with nonlinear kinematic and isotropic hardening and application to sheet metal forming. Int J Plast 26(5):659-687. https://doi.org/10.1016/j.ijplas.2009.09.008

24. Clifton RJ (1972) On the Equivalence of FeFp and FpFe. J Appl Mech 39(1):287-289. https://doi.org/10.1115/1.3422634

25. Lubarda VA (1999) Duality in constitutive formulation of finite-strain elastoplasticity based on F=FeFp and F=FpFe decompositions. Int J Plast 15(12):1277-1290. https://doi.org/10.1016/S0749-6419(99)00039-X

26. Lubarda VA (1991) Constitutive analysis of large elasto-plastic deformation based on the multiplicative decomposition of deformation gradient. Int J Solids Struct 27(7):885-895. https://doi.org/10.1016/0020-7683(91)90022-8

27. Bertram A (1999) An alternative approach to finite plasticity based on material isomorphisms. Int J Plast 15(3):353-374. https://doi.org/10.1016/S0749-6419(98)00074-6

28. Bouby C, Kondo D, Saxcé G (2015) A comparative analysis of two formulations for non linear hardening plasticity models: application to shakedown analysis. Eur J Mech A Solids 53:48-61. https://doi.org/10.1016/j.euromechsol.2015.03.001

29. Armstrong PJ, Frederick CO (1966) A mathematical representation of the multiaxial Bauschinger effect. C.E.G.B. Report RD/B/N 731, Berkeley Nuclear Laboratories, Berkeley, UK

30. Lemaitre J, Chaboche J-L (1994) Mechanics of solid materials. Cambridge University Press, Cambridge

31. Svendsen B, Arndt S, Klingbeil D, Sievert R (1998) Hyperelastic models for elastoplasticity with non-linear isotropic and kinematic hardening at large deformation. Int J Solids Struct 35(25):3363-3389. https://doi.org/10.1016/S0020-7683(98)00009-2

32. Dettmer W, Reese S (2004) On the theoretical and numerical modelling of Armstrong-Frederick kinematic hardening in the finite strain regime. Comput Methods Appl Mech Eng 193(1):87-116. https://doi.org/10.1016/j.cma.2003.09.005

33. Vladimirov IN, Pietryga MP, Reese S (2008) On the modelling of non-linear kinematic hardening at finite strains with application to springback-Comparison of time integration algorithms. Int J Numer Meth Eng 75(1):1-28. https://doi.org/10.1002/nme.2234

34. Dogui A, Sidoroff F (1985) Kinematic hardening in large elastoplastic strain. Eng Fract Mech 21(4):685-695. https://doi.org/10.1016/0013-7944(85)90078-5

35. Tsakmakis C (1996) Kinematic hardening rules in finite plasticity Part I: a constitutive approach. Continuum Mech Thermodyn 8(4):215-231. https://doi.org/10.1007/s001610050040

36. Tsakmakis C (1996) Kinematic hardening rules in finite plasticity Part II: some examples. Continuum Mech Thermodyn 8(4):233-246. https://doi.org/10.1007/s001610050041

37. Tsakmakis C, Willuweit A (2004) A comparative study of kinematic hardening rules at finite deformations. Int J Non-Linear Mech 39(4):539-554. https://doi.org/10.1016/S0020-7462(02)00221-4

38. Kratochvíl J (1973) On a finite strain theory of elastic-inelastic materials. Acta Mech 16(1):127-142. https://doi.org/10.1007/BF01177131

39. Gurtin ME, Anand L (2005) The decomposition F=FeFp, material symmetry, and plastic irrotationality for solids that are isotropic-viscoplastic or amorphous. Int J Plast 21(9):1686-1719. https://doi.org/10.1016/j.ijplas.2004.11.007

40. Loret B (1983) On the effects of plastic rotation in the finite deformation of anisotropic elastoplastic materials. Mech Mater 2(4):287-304. https://doi.org/10.1016/0167-6636(83)90021-2

41. Dafalias YF (1983) Corotational Rates for Kinematic Hardening at Large Plastic Deformations. J Appl Mech 50(3):561-565. https://doi.org/10.1115/1.3167091

42. Dafalias YF (1993) On multiple spins and texture development. Case study: kinematic and orthotropic hardening. Acta Mech 100(3):171-194. https://doi.org/10.1007/BF01174788

43. Dafalias YF (1984) The plastic spin concept and a simple illustration of its role in finite plastic transformations. Mech Mater 3(3):223-233. https://doi.org/10.1016/0167-6636(84)90021-8

44. Ulz MH (2011) A finite isoclinic elasto-plasticity model with orthotropic yield function and notion of plastic spin. Comput Methods Appl Mech Eng 200(21):1822-1832. https://doi.org/10.1016/j.cma.2011.01.020

45. Halphen B, Son Nguyen Q (1975) Sur les matériaux standard généralisés. J de Mécanique 14:39-63

46. Itskov M, Aksel N (2004) A constitutive model for orthotropic elasto-plasticity at large strains. Arch Appl Mech 74(1):75-91. https://doi.org/10.1007/BF02637210

47. Hashiguchi K (2019) Multiplicative hyperelastic-based plasticity for finite elastoplastic deformation/sliding: a comprehensive review. Arch Comput Methods Eng 26(3):597-637. https://doi.org/10.1007/s11831-018-9256-5

48. Foerch R, Besson J, Cailletaud G, Pilvin P (1997) Polymorphic constitutive equations in finite element codes. Comput Methods Appl Mech Eng 141(3):355-372. https://doi.org/10.1016/S0045-7825(96)01111-5

49. Holthusen H, Rothkranz C, Lamm L, Brepols T, Reese S (2023) Inelastic material formulations based on a co-rotated intermediate configuration-application to bioengineered tissues. J Mech Phys Solids 172:105174. https://doi.org/10.1016/j.jmps.2022.105174

50. Zheng Q-S (1994) Theory of representations for tensor functions–a unified invariant approach to constitutive equations. Appl Mech Rev 47(11):545-587. https://doi.org/10.1115/1.3111066

51. Raoult A (1986) Non-polyconvexity of the stored energy function of a Saint Venant–Kirchhoff material. Aplikace Matematiky 31(6):417-419

52. Holzapfel GA (2000) Nonlinear solid mechanics: a continuum approach for engineering. Wiley, Chichester

53. Ciarlet PG (1988) Three-dimensional elasticity. Elsevier, Amsterdam

54. Wallin M, Ristinmaa M (2005) Deformation gradient based kinematic hardening model. Int J Plast 21(10):2025-2050. https://doi.org/10.1016/j.ijplas.2005.01.007

55. Vujšević L, Lubarda VA (2002) Finite-strain thermoelasticity based on multiplicative decomposition of deformation gradient. Theoret Appl Mech 28-29:379-399. https://doi.org/10.2298/TAM0229379V

56. Bertram A (2003) Finite thermoplasticity based on isomorphisms. Int J Plast 19(11):2027-2050. https://doi.org/10.1016/S0749-6419(03)00057-3

57. Lu SCH, Pister KS (1975) Decomposition of deformation and representation of the free energy function for isotropic thermoelastic solids. Int J Solids Struct 11(7):927-934. https://doi.org/10.1016/0020-7683(75)90015-3
58. Teodosiu C, Sidoroff F (1976) A theory of finite elastoviscoplasticity of single crystals. Int J Eng Sci 14(2):165-176. https://doi.org/10.1016/0020-7225(76)90085-9

59. Lubarda VA (2004) Constitutive theories based on the multiplicative decomposition of deformation gradient: thermoelasticity, elastoplasticity, and biomechanics. Appl Mech Rev 57(2):95-108. https://doi.org/10.1115/1.1591000

60. Boyce MC, Montagut EL, Argon AS (1992) The effects of thermomechanical coupling on the cold drawing process of glassy polymers. Polymer Engineering & Science 32(16):1073-1085. https://doi.org/10.1002/pen.760321605

61. Kamlah M, Tsakmakis C (1999) Use of isotropic thermoelasticity laws in finite deformation viscoplasticity models. Continuum Mech Thermodyn 11(2):73-88. https://doi.org/10.1007/s001610050104

62. Yu J-S, Maniatty AM, Knorr DB (1997) Model for predicting thermal stresses in thin polycrystalline films. J Mech Phys Solids 45(4):511-534. https://doi.org/10.1016/S0022-5096(96)00104-4

63. Bertram A, Krawietz A (2012) On the introduction of thermoplasticity. Acta Mech 223(10):2257-2268. https://doi.org/10.1007/s00707-012-0700-6

64. Cailletaud G, Saï K (1995) Study of plastic/viscoplastic models with various inelastic mechanisms. Int J Plast 11(8):991-1005. https://doi.org/10.1016/S0749-6419(95)00040-2

65. Cailletaud G (2017) Multi-mechanism modeling of inelastic material behavior. ISTE Ltd/John Wiley and Sons Inc, Hoboken, NJ

66. Frederick CO, Armstrong PJ (2007) A mathematical representation of the multiaxial Bauschinger effect. Mater High Temp 24:1-26. https://doi.org/10.1179/096034007X207589

67. Chaboche JL (2008) A review of some plasticity and viscoplasticity constitutive theories. Int J Plast 24:1642-1693. https://doi.org/10.1016/j.ijplas.2008.03.009

68. Besson J, Cailletaud G, Chaboche J-L, Forest S (2009) Non-linear mechanics of materials. Springer, Berlin Heidelberg

69. Saxcé G (1992) Une généralisation de l'inégalité de Fenchel et ses applications aux lois constitutives. C.R. Acad. Sci. Paris, vol 314

70. Giessen E (1991) Micromechanical and thermodynamic aspects of the plastic spin. Int J Plast 7(5):365-386. https://doi.org/10.1016/0749-6419(91)90010-V

71. Aravas N (1994) Finite-strain anisotropic plasticity and the plastic spin. Modell Simul Mater Sci Eng 2(3A):483-504. https://doi.org/10.1088/0965-0393/2/3A/005

72. Fares N, Dvorak GJ (1991) Large elastic-plastic deformations of fibrous metal matrix composites. J Mech Phys Solids 39(6):725-744. https://doi.org/10.1016/0022-5096(91)90022-G

73. Böhlke T, Bertram A, Krempl E (2003) Modeling of deformation induced anisotropy in free-end torsion. Int J Plast 19:1867-1884. https://doi.org/10.1016/S0749-6419(03)00043-3

74. Neto EAdS, Peric D, Owen DRJ (2008) Computational methods for plasticity: theory and applications. John Wiley & Sons, Chichester

75. Weber G, Anand L (1990) Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. Comput Methods Appl Mech Eng 79(2):173-202. https://doi.org/10.1016/0045-7825(90)90131-5

76. Miehe C (1996) Multisurface thermoplasticity for single crystals at large strains in terms of Eulerian vector updates. Int J Solids Struct 33(20):3103-3130. https://doi.org/10.1016/0020-7683(95)00274-X

77. Truesdell C, Noll W (1965) The non-linear field theories of mechanics. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-46015-9_1

78. Méric L, Poubanne P, Cailletaud G (1991) Single crystal modeling for structural calculations. Part 1: Model presentation. J Eng Mater Technol 113:162-170

79. Duchene L, Lelotte T, Flores P, Bouvier S, Habraken A (2008) Rotation of axes for anisotropic metal in FEM simulations. Int J Plast 24(3):397-427. https://doi.org/10.1016/j.ijplas.2007.03.015

80. Peeters B, Hoferlin E, Van Houtte P, Aernoudt E (2001) Assessment of crystal plasticity based calculation of the lattice spin of polycrystalline metals for FE implementation. Int J Plast 17(6):819-836. https://doi.org/10.1016/S0749-6419(00)00070-X

81. Le LT, Ammar K, Forest S (2020) Efficient simulation of single and poly-crystal plasticity based on the pencil glide mechanism. Comptes Rendus Mécanique 348:846-876. https://doi.org/10.5802/crmeca.44

82. Méric L, Cailletaud G (1991) Single Crystal Modeling for Structural Calculations: Part 2-Finite Element Implementation. J Eng Mater Technol 113(1):171-182. https://doi.org/10.1115/1.2903375

83. Schröder J, Wick T, Reese S, Wriggers P, Müller R, Kollmannsberger S, Kästner M, Schwarz A, Igelbüscher M, Viebahn N, Bayat HR, Wulfinghoff S, Mang K, Rank E, Bog T, D'Angella D, Elhaddad M, Hennig P, Düster A, Garhuom W, Hubrich S, Walloth M, Wollner W, Kuhn C, Heister T (2021) A selection of benchmark problems in solid mechanics and applied mathematics. Arch Comput Methods Eng 28:713-751. https://doi.org/10.1007/s11831-020-09477-3

84. Schröder J, Wriggers P, Kraus A, Viebahn N (2022) In: Schröder, J., Wriggers, P. (eds.) Novel Finite Elements - Mixed, Hybrid and Virtual Element Formulations at Finite Strains for 3D Applications, pp. 37-67. Springer, Cham. https://doi.org/10.1007/978-3-030-92672-4_2

85. Lagioia R, Panteghini A (2016) On the existence of a unique class of yield and failure criteria comprising Tresca, von Mises, Drucker-Prager, Mohr-Coulomb, Galileo-Rankine, Matsuoka-Nakai and Lade-Duncan. Proc R Soc A 472:20150713. https://doi.org/10.1098/rspa.2015.0713

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.