A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate

# A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate

Marcos Latorre Francisco J. Montáns Escuela Técnica Superior de Ingeniería Aeronáutica y del Espacio
Universidad Politécnica de Madrid
Plaza Cardenal Cisneros, 3, 28040-Madrid, Spain
###### Abstract

In this paper we present a new general framework for anisotropic elastoplasticity at large strains. The new framework presents the following characteristics: (1) It is valid for non-moderate large strains, (2) it is valid for both elastic and plastic anisotropy, (3) its description in rate form is parallel to that of the infinitesimal formulation, (4) it is compatible with the multiplicative decomposition, (5) results in a similar framework in any stress-strain work-conjugate pair, (6) it is consistent with the principle of maximum plastic dissipation and (7) does not impose any restriction on the plastic spin, which must be given as an independent constitutive equation. Furthermore, when formulated in terms of logarithmic strains in the intermediate configuration: (8) it may be easily integrated using a classical backward-Euler rule resulting in an additive update. All these properties are obtained simply considering a plastic evolution in terms of a corrector rate of the proper elastic strain. This formulation presents a natural framework for elastoplasticity of both metals and soft materials and solves the so-called rate issue.

###### keywords:
Anisotropic material, Elastic-plastic material, Finite strains, Equations, Plastic flow rule.
journal: International Journal of Plasticity

## 1 Introduction

Constitutive models and integration algorithms for infinitesimal elastoplasticity are well established Simo98 ; KojicBatheBook ; LublinerBook . The currently favoured algorithmic formulations, either Cutting Plane Algorithms or Closest Point Projection ones are based on the concept of trial elastic predictor and subsequent plastic correction SimoHughesBook . The implementations of the most efficient closest point projection algorithms perform both phases in just two subsequent substeps MinCamMon . From the 70’s, quite a high number of formulations have been proposed to extend both the continuum and the computational small strain formulations to the finite deformation regime. Very different ingredients have been employed in these formulations, as for example different kinematic treatments of the constitutive equations, different forms of the internal elastic-plastic kinematic decomposition, different types of stress and strain measures being used, different internal variables chosen as the basic ones and, most controversially, different evolution equations for the plastic flow. The combinations of these ingredients have resulted into very different extended formulations ShutovIhlemann14 . However, as a common characteristic, all the formulations are developed with the main aim of preserving as much as possible the simplicity of the classical return mapping schemes of the infinitesimal theory Wilkins64 ; MaenchenSacks64 ; KriegKey76 through an algorithm that computes the closest point projection of the trial stresses onto the elastic domain.

The first strategies to model finite strain elastoplasticity were based on both an additive decomposition of the deformation rate tensor into elastic and plastic contributions and a hypoelastic relation for stresses TruesdellNollBook , see for example HibbittEtAl70 ; McMeekingRice75 ; KeyKrieg82 ; TaylorBecker83 among many others. Since the elastic stress relations are directly given in rate form and do not derive in general from a stored energy potential, some well-known problems may arise in these rate-form formulations, e.g. lack of objectivity of the resulting integration algorithms and the appearance of nonphysical energy dissipation in closed elastic cycles SimoPister84 ; KojicBathe87 . Incrementally objective integration algorithms HughesWinget80 ; PinskyOrtizPister83 overcome the former drawback; the selection of the proper objective stress rate, i.e. the corotational logarithmic rate in the so-called self-consistent Eulerian model XiaoBruhnsMeyers97 ; BruhnsXiaoMeyers99 ; XiaoBruhnsMeyers00 circumvents the latter one BrepolsVladimirovReese14 . Even though this approach is still being followed by several authors Teeriaho13 ; XiaoEtAl16 ; ZhuEtAl16 and may still be found in commercial finite element codes, the inherent difficulty associated to the preservation of objectivity in incremental algorithms makes these models less appealing from a computational standpoint RubinsteinAtluri83 ; BrepolsVladimirovReese14 .

Shortly afterwards the intrinsic problems of hypoelastic rate models arose, several hyperelastic frameworks formulated relative to different configurations emerged ArgyrisDoltsinis79 ; SimoOrtiz85 . Green-elastic, non-dissipative stresses are derived in these cases from a stored energy function, hence elastic cycles become path-independent and yield no dissipation GabrielBathe95 . Furthermore, objectivity requirements are automatically satisfied by construction of the hyperelastic constitutive relations SimoHughesBook .

In hyperelastic-based models, the argument of the stored energy potential from which the stresses locally derive is an internal elastic strain variable that has to be previously defined from the total deformation. Two approaches are common when large strains are considered. On the one hand, metric plasticity models propose an additive split of a given Lagrangian strain tensor into plastic and elastic contributions GreenNaghdi65 . On the other hand, multiplicative plasticity models are based on the multiplicative decomposition of the total deformation gradient into plastic and elastic parts Lee69 . The main advantage of the former type is that the proposed split is parallel to the infinitesimal one, where the additive decomposition of the total strain into plastic and elastic counterparts is properly performed, so these models somehow retain the desired simplicity of the small strain plasticity models SimoOrtiz85 ; Miehe98 ; PapadopoulusLu98 . Another immediate consequence is that these models are readily extended in order to include anisotropic elasticity and/or plasticity effects PapadopoulusLu01 ; Miehe02 ; LobleinSchroderGruttmann03 ; SansourWagner03 ; Ulz09 . However, it is well known that add hoc decompositions in terms of plastic metrics do not represent correctly the elastic part of the deformation under general, non-coaxial elastoplastic deformations LublinerBook ; Miehe02 ; GreenNaghdi71 ; Schmidt05 , hence its direct inclusion in the stored energy function may be questioned. For example, it has been found that these formulations do not yield a constant stress response when a perfectly plastic isotropic material is subjected to simple shear, a behavior which may be questionable Itskov04 . Furthermore, it has been recently shown NeffGhiba16 that these formulations may even modify the ellipticity properties of the stored energy function at some plastic deformation levels, giving unstable elastic spring back computations as a result, which seems an unrealistic response. On the contrary, multiplicative plasticity models are micromechanically motivated from single crystal metal plasticity Taylor38 ; Rice71 . The elastic part of the deformation gradient accounts for the elastic lattice deformation and the corresponding strain energy may be considered well defined. As a result, the mentioned plastic shear and elastic spring back degenerate responses do not occur in these physically sound models Itskov04 ; NeffGhiba16 .

Restricting now our attention to the widely accepted hyperelasto-plasticity formulations based on the multiplicative decomposition of the deformation gradient Lee69 ; Mandel74 , further kinematic and constitutive modelling aspects have to be defined. On one side, even though spatial quadratic strain measures were firstly employed Simo88 , they proved not to be natural in order to preserve plastic incompressibility, which had to be explicitly enforced in the update of the intermediate placement SimoMiehe92 . The fact that logarithmic strain measures inherit some properties from the infinitesimal ones, e.g. additiveness (only within principal directions), material-spatial metric preservation, same deviatoric-volumetric projections, etc., along with the excellent predictions that the logarithmic strain energy with constant coefficients provided for moderate elastic stretches Anand79 ; Anand86 , see also Rice75 ; RolphBathe84 , motivated the consideration of the quadratic Hencky strain energy in isotropic elastoplasticity formulations incorporating either isotropic or combined isotropic-kinematic hardening WeberAnand90 ; EterovicBathe90 ; PericOwenHonnor92 ; CuitinoOrtiz92 ; Simo92 . Exact preservation of plastic volume for pressure insensitive yield criteria is readily accomplished in this case. Moreover, the incremental schemes written in terms of logarithmic strains preserve the desired structure of the standard return mapping algorithms of classical plasticity models Simo92 , hence providing the simplest computational framework suitable for geometrically nonlinear finite element calculations.

On the other side, even though the use of logarithmic strain measures in actual finite strain computational elastoplasticity models has achieved a degree of common acceptance, a very controversial aspect of the theory still remains. This issue is the specific form that the evolution equations for the internal variables should adopt and how they must be further integrated Bruhns15 , a topic coined as the “rate issue” by Simó Simo92 . This issue originates, indeed, the key differences between the existing models. In this respect, the selection of the basic internal variable, whether elastic or plastic, in which the evolution equation is written becomes fundamental in a large deformation context. Evidently, this debate is irrelevant in the infinitesimal framework, where both the strains and the strain rates are fully additive. Early works Eckart48 ; Besseling66 ; Leonov76 suggest that the same strain variable on which the material response depends, i.e. the internal elastic strains, should govern the internal dissipation Rubin16 . This argument seems also reasonable from a numerical viewpoint taking into account that in classical integration algorithms Wilkins64 ; MaenchenSacks64 ; KriegKey76 the trial stresses, which are elastic in nature and directly computed from the trial elastic strains, govern the dissipative return onto the elastic domain during the plastic correction substep. Following this approach, Simó Simo92 used a continuum evolution equation for associative plastic flow explicitly expressed in terms of the Lie derivative of the elastic left Cauchy–Green deformation tensor (taken as the basic internal deformation variable SimoMiehe92 ). He then derived an exponential return mapping scheme to yield a Closest Point Projection algorithm formulated in elastic logarithmic strain space identical in structure to the infinitesimal one, hence solving the “rate issue” Simo92 . However, the computational model is formulated in principal directions and restricted to isotropy, so arguably that debated issue was only partially solved. Extensions of this approach to anisotropy are scarce, often involving important modifications regarding the standard return mapping algorithms (cf. Rubin16 and references therein).

Instead, the probably most common approach when modeling large strain multiplicative plasticity in the finite element context lies in the integration of evolution equations for the plastic deformation gradient, as done originally by Eterovic and Bathe EterovicBathe90 and Weber and Anand WeberAnand90 . The integration is performed through an exponential approximation to the incremental flow rule Simo98 , so these formulations are restricted to moderately large elastic strains EterovicBathe90 ; CamineroMontansBathe11 , which is certainly a minor issue in metal plasticity. However we note that it may be relevant from a computational standpoint if large steps are involved because the trial substep may involve non-moderate large strains. Unlike Simó’s approach, these models retain a full tensorial formulation, so further consideration of elastic and/or plastic anisotropy is amenable ChattiEtAl01 ; HanEtAL03 ; EidelGruttman03 ; MenzelEtAl05 ; SansourKarsajSoric06 ; MontansBathe07 ; KimMontansBathe09 ; VladimirovPietrygaReese10 ; CamineroMontansBathe11 . However, the consideration of elastic anisotropy in these models has several implications in both the continuum and the algorithmic formulations, all of them derived from the fact that the resulting thermodynamical stress tensor in the intermediate configuration, i.e. the Mandel stress tensor Mandel72 , is non-symmetric in general. Interestingly, the symmetric part of this stress tensor is, in practice, work-conjugate of the elastic logarithmic strain tensor for moderately large elastic deformations, which greatly simplifies the algorithmic treatment CamineroMontansBathe11 in anisotropic metal plasticity applications. As a result, the model in CamineroMontansBathe11 , formulated in terms of generalized Kirchhoff stresses instead of Kirchhoff stresses and with the additional assumption of vanishing plastic spin, becomes the natural generalization of the Eterovic and Bathe model EterovicBathe90 to the fully anisotropic case, retaining at the same time the interesting features of the small strain elastoplasticity theory and algorithms.

Summarizing, the computational model of Caminero et al. CamineroMontansBathe11 is adequate for anisotropic elastoplasticity but not for non-moderate large elastic deformations. In contrast, the Simó formulation Simo92 is valid for large elastic strains but not for phenomenological anisotropic elastoplasticity. In this work we present a novel continuum elastoplasticity framework in full space description valid for anisotropic elastoplasticity and large elastic deformations consistent with the Lee multiplicative decomposition. The main novelty is that, generalizing Simó’s approach Simo92 , internal elastic deformation variables are taken as the basic variables that govern the local dissipation process. The dissipation inequality is reinterpreted taking into account that the chosen internal elastic tensorial variable depends on the respective internal plastic variable and also on the external one. In this reinterpretation we take special advantage of the concepts of partial differentiation and mapping tensors LatMonAPM2016 . The procedure is general and may be described in different configurations and in terms of different stress and strain measures, yielding as a result dissipation inequalities that are fully equivalent to each other. Respective thermodynamical symmetric stress tensors and associative flow rules expressed in terms of corrector elastic strain rates and general yield functions are trivially obtained consistently with the principle of maximum dissipation. We recover the Simó framework from our spatial formulation specialized to isotropy and with the additional assumption of vanishing plastic spin, as implicitly assumed in Ref. Simo92 , see also BonetWoodBook . Exactly as it occurs in the infinitesimal theory, in all the descriptions being addressed the plastic spin does not take explicit part in the associative six-dimensional flow rules being derived, hence bypassing the necessity of postulating a flow rule for the plastic spin as an additional hypothesis in the dissipation equation Lubliner86 . Special advantage is taken when the continuum formulation is written in terms of the logarithmic elastic strain tensor LatMonIJSS2014  and its work-conjugated generalized Kirchhoff symmetric stress tensor, both defined in the intermediate configuration. Then, the continuum formulation mimics the additive description in rate form of the infinitesimal elastoplasticity theory, the only differences coming from the additional geometrical nonlinearities arising in a finite deformation context. Furthermore, the unconventional appearance Simo92 of the well-known continuum evolution equation defining plastic flow in terms of the Lie derivative of the elastic left Cauchy–Green tensor in the current configuration BonetWoodBook makes way for a conventional evolution equation in terms of the elastic logarithmic strain rate tensor in the intermediate placement, hence simplifying the continuum formulation to a great extent and definitively solving the “rate issue” directly in the logarithmic strain space. Remarkably, with the present multiplicative elastoplasticity model at hand, the generally non-symmetric stress tensor that has traditionally governed the plastic dissipation in the intermediate configuration, i.e. the Mandel stress tensor, is no longer needed.

The rate formulation that we present herein in terms of logarithmic strains in the intermediate configuration may be immediately recast in a remarkably simple incremental form by direct backward-Euler integration which results in integration algorithms of similar additive structure to those of the infinitesimal framework. Indeed, the formulation derived herein is equivalent in many aspects to the anisotropic finite strain viscoelasticity model based on logarithmic strains and the Sidoroff multiplicative decomposition that we presented in Ref. LatMonCM2015 . As done therein, a first order accurate backward-Euler algorithm could be directly employed over the corrector logarithmic elastic strain rate flow rule obtained herein to yield a return mapping scheme in full tensorial form, valid for anisotropic finite strain responses, that would preserve the appealing structure of the classical return mapping schemes of infinitesimal plasticity without modification. For the matter of simplicity in the exposition of the new elastoplasticity framework, we do not include kinematic hardening effects in the formulation. Nevertheless, its further consideration would be straightforward.

The rest of the paper is organized as follows. We next present in Section 2 the ideas for infinitesimal elastoplasticity in order to motivate and to prepare the parallelism with the finite strain formulation. Thereafter we present in Section 3 the large strain formulation in the spatial configuration performing such parallelism. We then particularize the present proposal to isotropy and demonstrate that some well-known formulations which are restricted to isotropy are recovered as a particular case from the more general, but at the same time simpler, anisotropic one. Section 4 is devoted to the formulation in the intermediate configuration, where a comparison with existing formulations is presented and some difficulties encountered in the literature are discussed. Section 5 presents the new approach to the problem at the intermediate configuration, both for quadratic strain measures and for our favoured logarithmic ones. In that section we also discuss the advantages and possibilities of the present framework.

## 2 Infinitesimal elastoplasticity: two equivalent descriptions

The purpose of this section is to motivate the concepts in the simpler infinitesimal description, showing a new subtle view of these equations which, thereafter result in a remarkable parallelism with the large strain formulations.

Consider the Prandtl (friction-spring) rheological model for small strains shown in Figure 1

where and are the external, measurable infinitesimal strains and engineering stresses, respectively, and and are internal, non-measurable infinitesimal strains describing the internal elastic and plastic behaviors. The internal strains relate to the external ones through

 \boldmathε \unboldmath=\boldmathε % \unboldmathe+\boldmathε \unboldmathp (1)

so if we know the total deformation and one internal variable, then the other internal variable is uniquely determined. We will consider and as the independent variables of the dissipative system and will be the dependent internal variable. The following two-variable dependence emerges for

 \boldmathε \unboldmathe(\boldmathε \unboldmath,\boldmathε \unboldmathp)=\boldmathε \unboldmath−\boldmathε \unboldmathp (2)

which provides also a relation between the corresponding strain rate tensors—we use the notation for partial differentiation

 \boldmath˙ε \unboldmathe=∂\boldmathε \unboldmathe∂\boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath:\boldmath˙ε \unboldmath+∂\boldmathε \unboldmathe∂\boldmathε \unboldmathp∣∣∣\boldmath˙ε \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath:\boldmath˙ε \unboldmathp=IS:\boldmath˙ε \unboldmath−IS:% \boldmath˙ε \unboldmathp=\boldmath˙ε \unboldmath−\boldmath˙ε \unboldmathp=\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmathp=% \boldmath0 \unboldmath\boldmath \unboldmath+\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmath=\boldmath% 0 \unboldmath\boldmath \unboldmath (3)

where stands for the fourth-order (symmetric) identity tensor

 (IS)ijkl=12(δikδjl+δilδjk) (4)

For further use, we define the following partial contributions to the elastic strain rate tensor

 \boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmathp=% \boldmath0 \unboldmath\boldmath \unboldmath=∂\boldmathε \unboldmathe∂\boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 % \unboldmath\boldmath \unboldmath:\boldmath% ˙ε \unboldmath=IS:\boldmath˙ε \unboldmath=\boldmath˙ε \unboldmath (5)

and

 \boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmath=\boldmath% 0 \unboldmath\boldmath \unboldmath=∂\boldmathε \unboldmathe∂% \boldmathε \unboldmathp∣∣∣% \boldmath˙ε \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath:\boldmath˙ε \unboldmathp=−IS:\boldmath˙ε \unboldmathp=−\boldmath˙ε % \unboldmathp (6)

The stored energy in the device of Figure 1 is given in terms of the internal elastic deformation, i.e. . The (non-negative) dissipation rate is calculated from the stress power and the total strain energy rate through

 D=P−˙Ψ≥0 (7)

which can be written as

 D=\boldmathσ \unboldmath:\boldmath˙ε \unboldmath−\boldmathσ \unboldmath|e:\boldmath˙ε \unboldmathe≥0 (8)

where we have introduced the following notation for the total gradient—we use the notation for total differentiation

 \boldmathσ \unboldmath|e:=dΨ(% \boldmathε \unboldmathe)d\boldmathε \unboldmathe (9)

No dissipation takes place if we consider an isolated evolution of the external, independent variable without internal variable evolution, i.e. with . Then, from Eq. (3) and Eq. (8) reads

 D=\boldmathσ \unboldmath:\boldmath˙ε \unboldmath−\boldmathσ \unboldmath|e:\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmathp=% \boldmath0 \unboldmath\boldmath \unboldmath=(\boldmathσ \unboldmath−\boldmathσ % \unboldmath|e:∂\boldmathε % \unboldmathe∂\boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath):\boldmath˙ε \unboldmath=0\quad if% \boldmath˙ε \unboldmathp=\boldmath0 \unboldmath (10)

which yields

 \boldmathσ \unboldmath=\boldmathσ \unboldmath|e:∂\boldmathε \unboldmathe∂\boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 % \unboldmath\boldmath \unboldmath=\boldmath% σ \unboldmath|e:IS=\boldmathσ % \unboldmath|e (11)

where we recognize the following definition based on a chain rule operation—note the abuse of notation ; we keep the dependencies explicitly stated when the distinction is needed

 σ =\boldmathσ \unboldmath|e:∂\boldmathε \unboldmathe∂% \boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath=dΨd\boldmathε \unboldmathe:∂\boldmathε \unboldmathe∂\boldmathε % \unboldmath∣∣∣\boldmath˙ε % \unboldmathp=\boldmath0 \unboldmath% \boldmath \unboldmath =dΨ(\boldmathε \unboldmathe)d\boldmathε \unboldmathe:∂\boldmathε \unboldmathe(\boldmathε \unboldmath,\boldmathε \unboldmathp)∂\boldmathε \unboldmath=∂Ψ(\boldmathε \unboldmath,% \boldmathε \unboldmathp)∂\boldmathε \unboldmath≡∂Ψ∂% \boldmathε \unboldmath∣∣∣\boldmath˙ε \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath (12)

These definitions based on the concept of partial differentiation relate internal variables with external ones from a purely kinematical standpoint and will prove extremely useful in the finite deformation context, where they will furnish the proper pull-back and push-forward operations between the different configurations being defined.

Consider now an isolated variation of the other independent variable in the problem, i.e. the case for which and , which note is a purely internal (dissipative) evolution. Then from Eq. (3) . The dissipation inequality of Eq. (8) must be positive because plastic deformation is taking place

 D=−\boldmathσ \unboldmath|e:% \boldmath˙ε \unboldmathe∣∣% \boldmath˙ε \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath>0\quad if% \boldmath˙ε \unboldmathp≠\boldmath0 % \unboldmath (13)

We arrive at the same expression of Eq. (13) if we consider the most general case for which both independent variables are simultaneously evolving, i.e. and . Hence note that both Eqs. (10) and (13) hold if either or , so only the respective condition over is indicated in those equations. Since in the infinitesimal framework of this section and , recall Eqs. (11) and (6), just in this case we can write Eq. (13) in its conventional form

 D=\boldmathσ \unboldmath:\boldmath˙ε \unboldmathp>0\quad if\boldmath˙ε \unboldmathp≠\boldmath0 \unboldmath (14)

i.e. the dissipation must be positive when the (six-dimensional) frictional element in Figure 1 experiences slip. Interestingly, Equations (13) and (14) represent both the same physical concept, the former written in terms of the partial contribution to the rate of the dependent internal variable and the latter written in terms of the total rate of the independent internal variable . However note that they present a clearly different interpretation which will become relevant in the large strain framework.

### 2.1 Local evolution equation in terms of \boldmath˙ε \unboldmathe|% \boldmath˙ε \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath

Equation (13) is automatically fulfilled if we choose the following evolution equation for the internal strains

 −\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmath=\boldmath% 0 \unboldmath\boldmath \unboldmath=˙γ1kN:\boldmathσ \unboldmath|e (15)

which yields

 D=\boldmathσ \unboldmath|e:N:\boldmathσ \unboldmath|ek2k˙γ>0% \quad if\boldmath˙ε \unboldmathp≠0 (16)

where is a fully symmetric positive definite fourth-order tensor, is the characteristic yield stress of the internal frictional element of Figure 1 and is the plastic strain rate component which is power-conjugate of the stress-like variable , as we see just below. If the internal yield stress is constant, the model describes the perfect plasticity case. If increases with an increment of the amount of plastic deformation, namely , the model may incorporate non-linear isotropic hardening effects. We rephrase the dissipation Equation (16) as

 D=1k2(\boldmathσ \unboldmath|e:N:\boldmathσ \unboldmath|e−k2)k˙γ+k˙γ>0\quad if˙γ>0 (17)

then we immediately recognize the yield function and the loading-unloading conditions

 (18)

and

 f(\boldmathσ \unboldmath|e,k)=\boldmathσ % \unboldmath|e:N:\boldmathσ \unboldmath|e−k2<0\quad⇒˙γ=0 (19)

so we obtain the plastic dissipation (if any) as given by the (scalar) flow stress times the (scalar) frictional strain rate for .

Equation (15) may be reinterpreted in terms of the yield function gradient to give the following associative flow rule for the internal elastic strains evolution

 −\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmath=\boldmath% 0 \unboldmath\boldmath \unboldmath=˙γ1k∇ϕ (20)

where we have introduced the quadratic form

 ϕ(\boldmathσ \unboldmath|e)=12\boldmath% σ \unboldmath|e:N:\boldmathσ \unboldmath|e (21)

for the matter of notation simplicity, so and .

### 2.2 Local evolution equation in terms of \boldmath˙ε \unboldmathp

Using the equivalences given in Eqs. (11) and (6), the yield function of Eq. (18) is given in terms of the (external) stress tensor as

 f(\boldmathσ \unboldmath,k)=\boldmathσ % \unboldmath:N:\boldmathσ \unboldmath−k2=2ϕ(\boldmathσ \unboldmath)−k2=0\quad if% ˙γ>0 (22)

and the associative flow rule of Eq. (20) adopts the usual expression in terms of the (internal) plastic strain rate tensor , cf. Eq. (2.5.6) of Ref. SimoHughesBook or Eq. (87) of Ref. MinCamMon

 \boldmath˙ε \unboldmathp=˙γ∇ϕ√\boldmathσ \unboldmath:N:% \boldmathσ \unboldmath (23)

As we discuss below, the interpretation given in Eq. (20) greatly facilitates the extension of the infinitesimal formulation to the finite strain context without modification.

### 2.3 Description in terms of trial and corrector elastic strain rates

It is apparent from the foregoing results that, in practice, no distinction is needed within the infinitesimal framework regarding both the selection of either or as the basic internal deformation variable and the selection of either or as the basic stress tensor. In what follows, however, we keep on developing the infinitesimal formulation in terms of and , which will let us take special advantage of the functional dependencies and .

Regarding the evolution of elastic variables, whether strains or stresses, it is convenient to introduce the concepts of trial and corrector elastic strain rates in Eq. (3). This decomposition in rate form is the origin of the trial elastic predictor, for which is frozen, and plastic corrector, for which is frozen, operator split typically employed for elastic internal variables in computational inelasticity within an algorithmic framework. Accordingly, we define within the continuum theory

 \boldmath˙ε \unboldmathe=\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath+\boldmath˙ε \unboldmathe∣∣\boldmath˙ε \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath=:tr\boldmath˙ε% \unboldmathe+ct\boldmath˙ε \unboldmathe (24)

where the superscripts tr and ct stand for trial and corrector respectively. Interestingly, the concepts of trial and corrector elastic rates emerge in the finite deformation multiplicative framework developed below without modification with respect to the infinitesimal case, so we will be able to directly compare the small and large strain formulations equation by equation. We note that elastoplasticity models based on plastic metrics have traditionally followed the same philosophy, but departing from the standard rate decomposition

 \boldmath˙ε \unboldmathe=\boldmath˙ε \unboldmath−\boldmath˙ε \unboldmathp (25)

which, however, leads to additive Lagrangian formulations GreenNaghdi65 , Miehe98 , PapadopoulusLu98 , PapadopoulusLu01 , Miehe02 , LobleinSchroderGruttmann03 , SansourWagner03 , Ulz09 that are not generally consistent with the finite strain multiplicative decomposition, as it is well-known GreenNaghdi71 , Itskov04 , Schmidt05 , NeffGhiba16 .

For further comparison, we rephrase both the dissipation inequality of Eq. (13) and the associative flow rule of Eq. (20) in terms of the corrector elastic strain rate as

 D=−\boldmathσ \unboldmath|e:ct% \boldmath˙ε \unboldmathe>0\quad if˙γ>0 (26)

and

 ct\boldmath˙ε \unboldmathe=−˙γ1k∇ϕ (27)

Note that the elastic strain correction performed in CPP algorithms and defined in Eq. (27) enforce the instantaneous closest point projection onto the elastic domain, i.e. the normality rule in the continuum setting.

In the case we do not consider a potential, then the formulation is usually referred to as generalized plasticity PastorZ , which is a generalization of nonassociative plasticity typically used in soils Borjabook . However, we can alternatively take

 ct\boldmath˙ε \unboldmathe=−˙γ1k\boldmathG \unboldmath(\boldmathσ % \unboldmath|e) (28)

where the prescribed second-order tensor function defines the direction of plastic flow. So Eq. (26) reads

 D=˙γ1k\boldmathσ \unboldmath|e:\boldmathG \unboldmathif˙γ>0 (29)

even though positive dissipation and a fully symmetric linearization of the continuum theory are not guaranteed in this case Simo98 . Note that for associative plasticity.

### 2.4 Maximum Plastic Dissipation

We assume now the existence of another arbitrary stress field different from the actual stress field , as given in Eq. (11). The dissipation originated by during the same plastic flow process would be—cf. Eq. (26)

 DΣ=−\boldmathΣ \unboldmath|e:ct\boldmath˙ε \unboldmathe\quad if˙γ>0 (30)

The evolution of plastic flow, e.g. Eq. (27), is said to obey the Principle of Maximum Dissipation if

 D−DΣ>0 (31)

for any admissible stress field , i.e. with . Considering the associative flow rule of Eq. (27), we arrive at

 D−DΣ=−(\boldmathσ \unboldmath|e−\boldmathΣ \unboldmath|e):ct\boldmath˙ε \unboldmathe=˙γ1k(\boldmathσ \unboldmath|e−\boldmathΣ \unboldmath|e):∇ϕ (32)

If is a strictly convex function and is admissible

 D−DΣ=˙γ1k(\boldmathσ \unboldmath|e−\boldmathΣ \unboldmath|e):∇ϕ=˙γ12k(\boldmathσ \unboldmath|e−\boldmathΣ \unboldmath|e):∇f>0 (33)

i.e. maximum dissipation in the system is guaranteed (the equal sign would be possible if non-strictly convex functions are considered, as for example Tresca’s one).

In all the finite strain cases addressed below if the corresponding associative flow rule for each case is considered. Indeed, this principle must hold in any arbitrary stress-strain work-conjugate couple, but if guaranteed in one of them, will hold in any of them by invariance of power.

## 3 Finite strain anisotropic elastoplasticity formulated in the current configuration

We present in this section a new framework for finite strain anisotropic elastoplasticity formulated in the current configuration in which the basic internal variables are elastic in nature. Once the corresponding dependencies are identified, the theory is further developed taking advantage of the previously introduced concepts of partial differentiation, mapping tensors and the trial-corrector decomposition of internal elastic variables in rate form. With the exception of the geometrical nonlinearities being introduced, the formulation yields identical expressions to those derived above for infinitesimal plasticity.

### 3.1 Multiplicative decomposition

The so-called Lee multiplicative decomposition Lee69  states the decomposition of the deformation gradient into an elastic part and a plastic part as

 \boldmathX \unboldmath=\boldmathX \unboldmathe\boldmathX \unboldmathp (34)

When using this decomposition, a superimposed rigid body motion by an orthogonal proper tensor results into

 \boldmathX \unboldmath+=\boldmathQX \unboldmath=\boldmathX \unboldmath+e\boldmathX \unboldmath+p=(\boldmathQX \unboldmathe)(% \boldmathX \unboldmathp) (35)

so the rigid body motion naturally enters the “elastic” gradient, whereas the plastic gradient remains unaltered. A much debated issue is the uniqueness of the intermediate configuration arising from since any arbitrary rotation tensor with its inverse may be inserted such that , so the decomposition of Eq. (34) is unique up to a rigid body rotation of the intermediate configuration. However, in practice, since is path dependent and is integrated step-by-step in an incremental fashion in computational elastoplasticity algorithms CamineroMontansBathe11 ; MontansBenitezCaminero12 , we consider that it is uniquely determined at all times.

### 3.2 Trial and corrector elastic deformation rate tensors

Consider the following additive decomposition of the spatial velocity gradient tensor

 \boldmathl \unboldmath:=\boldmath˙XX \unboldmath−1=\boldmath˙X \unboldmathe\boldmathX % \unboldmath−1e+\boldmathX \unboldmathe% \boldmath˙X \unboldmathp\boldmathX \unboldmath−1p\boldmathX \unboldmath−1e=\boldmathl % \unboldmathe+\boldmathX \unboldmathe\boldmathl% \unboldmathp\boldmathX \unboldmath−1e (36)

where we define the elastic and plastic velocity gradients as

 \boldmathl \unboldmathe:=\boldmath˙X \unboldmathe\boldmathX \unboldmath−1e\quad and\quad%\boldmath$l$\unboldmathp:=\boldmath˙X \unboldmathp\boldmathX \unboldmath−1p (37)

We note that lies in the spatial configuration, whereas operates in the intermediate configuration. The deformation rate tensor (the symmetric part of ) and the spin tensor (its skew-symmetric part) are

 \boldmathd \unboldmath=sym(\boldmathl \unboldmath)\quad and\quad\boldmathw \unboldmath=skw(\boldmathl \unboldmath) (38)

The elastic and plastic velocity gradient tensors also admit the corresponding decomposition into deformation rate and spin counterparts, and , thereby from Eq. (36)

 \boldmathd \unboldmath=\boldmathd \unboldmathe+sym(\boldmathX \unboldmathe\boldmathl \unboldmath% p\boldmathX \unboldmath−1e) (39)
 \boldmathw \unboldmath=\boldmathw \unboldmathe+skw(\boldmathX \unboldmathe\boldmathl \unboldmath% p\boldmathX \unboldmath−1e) (40)

In general, from Eq. (39) we can consider the elastic deformation rate tensor as a two-variable function of the deformation rate tensor and the plastic velocity gradient tensor (including the plastic spin ) through

 \boldmathd \unboldmathe(\boldmathd,l \unboldmathp)=\boldmathd \unboldmath−sym(\boldmathX % \unboldmathe\boldmathl \unboldmathp\boldmathX % \unboldmath−1e) (41)

which can be expressed in the following rate-form formats—compare with Eqs. (3) and (24)

 \boldmathd \unboldmathe=Mded∣∣\boldmathl \unboldmathp=\boldmath0 % \unboldmath\boldmath \unboldmath:\boldmath% d \unboldmath+Mdelp∣∣%\boldmath$d$\unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath:\boldmathl \unboldmathp=\boldmathd \unboldmathe∣∣% \boldmathl \unboldmathp=\boldmath0 \unboldmath\boldmath \unboldmath+\boldmathd % \unboldmathe∣∣\boldmathd \unboldmath=\boldmath0 \unboldmath\boldmath \unboldmath=tr\boldmathd \unboldmathe+ct\boldmathd % \unboldmathe (42)

where and are mapping tensors CamineroMontansBathe11 ; LatMonAPM2016 which allow us to define the following partial contributions to the elastic deformation rate tensor

 tr\boldmathd \unboldmathe:=Mded∣∣\boldmathl \unboldmathp=\boldmath0% \unboldmath\boldmath \unboldmath:% \boldmathd \unboldmath=IS:\boldmathd \unboldmath=\boldmathd \unboldmath (43)

and

 ct\boldmathd \unboldmathe=Mdelp∣∣\boldmathd \unboldmath=\boldmath0 % \unboldmath\boldmath \unboldmath:\boldmath% l \unboldmathp=−12(\boldmathX \unboldmathe⊙\boldmathX \unboldmath−Te+\boldmathX % \unboldmath−Te⊡\boldmathX \unboldmathe):\boldmathl \unboldmathp=−sym(\boldmathX % \unboldmathe\boldmathl \unboldmathp\boldmathX % \unboldmath−1e) (44)

with and .

It is frequently assumed in computational plasticity that the plastic spin vanishes, namely , so its effects in the dissipation inequality are not taken into account. However, as in the small strain case discussed above, the plastic spin evolves independently of the normality flow rules being developed below in terms of corrector elastic rates, so no additional assumptions over will be prescribed by the dissipation process Lubliner86 . The a priori undetermined intermediate configuration, defined by , would become determined once an independent constitutive equation for the plastic spin is specified Simo98 , MontansBathe07 , KimMontansBathe09 , which is strictly needed in order to complete the model formulation.

### 3.3 Dissipation inequality and flow rule in terms of ct\boldmathd \unboldmathe

From purely physical grounds, we know that the strain energy function locally depends on an elastic measure of the deformation. Hence it may be expressed in terms of a Lagrangian-like elastic strain tensor lying in the intermediate configuration, e.g. the elastic Green–Lagrange-like strains where is the second-order identity tensor, as

 ΨA=ΨA(\boldmathA \unboldmathe,\boldmath% a \unboldmath1⊗\boldmatha \unboldmath1,% \boldmatha \unboldmath2⊗\boldmatha \unboldmath2) (45)

where we have additionally assumed that the material is orthotropic, with and (and ) defining the orthogonal preferred directions in the intermediate configuration. As a first step in the derivation of more complex formulations including texture evolution, which involves an experimentally motivated constitutive equation additional to that for , see examples in Ref. KimMontansBathe09 and references therein, we assume in this work that the texture of the material is permanent and independent of the plastic spin. That is, we consider the case for which is given as an additional equation so that the Lee decomposition is completely defined at each instant and we take as a simplifying assumption for the stresses update. Subsequently, the material time derivative of the Lagrangian potential may be expressed in terms of variables lying in the current configuration through

 ˙ΨA=dΨA(\boldmathA \unboldmathe)d\boldmathA \unboldmathe:\boldmath˙A % \unboldmathe=\boldmathS \unboldmath|e:\boldmath˙A \unboldmathe=\boldmathS \unboldmath|e:% \boldmathX \unboldmathTe⊙\boldmathX \unboldmathTe:\boldmathd \unboldmathe=\boldmathτ % \unboldmath|e:\boldmathd \unboldmathe (46)

where we have used the purely kinematical pull-back operation over (lying in the current configuration) that gives (lying in the intermediate configuration) —see LatMonAPM2016

 \boldmath˙A \unboldmathe=\boldmathX \unboldmathTe\boldmathd \unboldmathe\boldmathX % \unboldmathe=\boldmathX \unboldmathTe⊙% \boldmathX \unboldmathTe:\boldmathd \unboldmathe=:M˙Aede:\boldmathd \unboldmathe (47)

which provides as a result the also purely kinematical push-forward operation over the internal elastic second Piola–Kirchhoff stress tensor (lying in the intermediate configuration)

 \boldmathS \unboldmath|e:=dΨA(\boldmath% A \unboldmathe)d\boldmathA \unboldmathe (48)

that gives the internal elastic Kirchhoff stress tensor (lying in the current configuration)

 \boldmathτ \unboldmath|e:=\boldmathS \unboldmath|e:M˙Aede