Plane-stress constrained multiplicative hyperelasto-plasticity with nonlinear kinematic hardening. Consistent theory based on elastic corrector rates and algorithmic implementation

Journal Pre-proof Plane-stress constrained multiplicative hyperelasto-plasticity with nonlinear kinematic hardening. Consistent theory based on elastic corrector rates and algorithmic implementation. K. Nguyen, Miguel A. Sanz, Francisco J. Montáns PII:

S0749-6419(19)30336-5

DOI:

https://doi.org/10.1016/j.ijplas.2019.08.017

Reference:

INTPLA 2592

To appear in:

International Journal of Plasticity

Received Date: 7 May 2019 Revised Date:

7 August 2019

Accepted Date: 23 August 2019

Please cite this article as: Nguyen ,, K., Sanz, M.A., Montáns, F.J., Plane-stress constrained multiplicative hyperelasto-plasticity with nonlinear kinematic hardening. Consistent theory based on elastic corrector rates and algorithmic implementation., International Journal of Plasticity, https:// doi.org/10.1016/j.ijplas.2019.08.017. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.

Plane-stress constrained multiplicative hyperelasto-plasticity with nonlinear kinematic hardening. Consistent theory based on elastic corrector rates and algorithmic implementation. K. Nguyena , Miguel A. Sanza , Francisco J. Mont´ansa,∗ a

Escuela T´ecnica Superior de Ingenier´ıa Aeron´ autica y del Espacio, Universidad Polit´ecnica de Madrid, Pza. Cardenal Cisneros, 28040, Madrid

Abstract Small-strain plane-stress projected formulations are typically employed in 2D, plate and thin membrane problems because of their efficiency. However, at large strains, equivalent algorithms have been precluded by the difficulties of kinematics, specially using the multiplicative decomposition and nonlinear kinematic hardening. In this paper, we present a plane-stress constrained formulation for large strain multiplicative hyperelasto-plasticity based on the novel framework which uses elastic strain corrector rates. The proposed plane-stress projected model is consistent with the principle of maximum dissipation as in the 3D/plane-strain case, and it is valid for any stored energy function (e.g. for soft materials). Hyperelasticity-based nonlinear kinematic hardening at large strains, which preserves Masing’s rules if desired, is also included in the formulation. The consistent tangent moduli tensor for the plane stress subspace is also derived to provide the asymptotic second order convergence of Newton algorithms during global equilibrium iterations. Some numerical examples are presented in order to evaluate both the congruency with the equivalent 3D formulation and the numerical performance of the proposed plane-stress projected algorithm. Keywords: large strains, plane stress plasticity, logarithmic strains, multiplicative decomposition, nonlinear kinematic hardening.

1. Introduction In the analysis of plates and thin membranes it is frequent to use the condition of plane stress, meaning that the out-of-plane stress components vanish (i.e. they are neglected when compared to the in-plane components). This assumption results in very important savings in memory and in computational time [1] so it is systematically employed in linear elasticity. In a finite element program, plane strain and axi-symmetric cases present no difficulty because since finite element procedures are strain-driven, they are considered in a straightforward ∗

Corresponding author Email addresses: [email protected] (K. Nguyen), [email protected] (Miguel A. Sanz), [email protected] (Francisco J. Mont´ ans )

Preprint submitted to Elsevier

September 6, 2019

manner by the simple consideration of the relevant (unmodified in-plane) equations from the 3D case. However, in the plane stress case, the plane stress constraint implies that the problems have mixed unknowns, namely the in-plane stresses and the out-of-plane strains. In linear elasticity, the plane stress conditions are easily imposed at the Gauss point level through a condensed matrix of elastic constants, so the rest of the finite element formulation is identical to the plane strain case [1]. Unfortunately, in the nonlinear case the situation is more complex. There are mainly three approaches to enforce the plane stress constraints in elastoplasticity, see reviews [2, 3]. One option is to impose the plane stress constraint (usually in the third dimension) in the 3D equations through a nested loop at the Gauss-point level, or alternatively through other technique like Lagrange multipliers. This approach has been used in [4, 5, 6, 7], among others, and has the advantage of being a simple modification of the 3D theory, but the disadvantage of being expensive [3]. The second option is to impose the constraints at the structural level, as proposed in [8]. The main advantage of this approach is that the material routines are unmodified respect to the 3D case, but this approach results also in little gain over a 3D formulation, so it is rarely used. The third option is a planestress projected formulation. In this latter option, all the formulation is developed already constrained to the plane-stress space at the Gauss-point level, and involves only the in-plane stresses and strains; out-of-plane strains may be computed aposteriori. This type of approach has an important popularity in the scientific community (see [9, 10, 11, 12, 13, 3, 14], among others), because it is a substantially more efficient procedure and it is in essence what it is done in plane stress linear elasticity. However, the development of this class of formulations is much more involved, so few algorithms of this type have been presented even for infinitesimal elastoplasticity, being the best known the algorithm of Simo and Taylor [15, 14]. Other algorithms are those of Jetteur [16] and Ramm and Matzenmiller [17]. For the case of nonlinear kinematic hardening of the Armstrong-Frederick type [18, 19], but in the infinitesimal strains context, Ohno et al [20] and Hopperstad and Remseth [21] have provided plane-stress projected algorithms. For multiplicative elastoplasticity, due to the inherent difficulties imposed by nonlinear kinematics (difficulties already present for developing sound 3D formulations), to the authors’ knowledge a plane-stress projected elastoplastic formulation is available only for isotropically hardened plasticity in principal logarithmic strains, and does not include nonlinear kinematic hardening nor arbitrary stored energy functions nor nonmoderate elastic strains, see e.g. [3] and therein references, and a similar approach in [22]. Hence, the purpose of the present paper is to fill this important gap in a consistent manner with multiplicative hyperelasto-plasticity, including the nonlinear kinematic hardening for nonproportional and cyclic loading cases. Efficient finite element simulations which use large strain continuum plasticity models are now systematically employed and required in analysing many processes like manufacturing, fatigue, ductile fracture, etc [23, 24, 25, 1]. Plane stress cases are abundant, for example, in sheet forming and crash-worthiness, where new metallic and polymeric materials are being employed. Then, large strain formulations for plasticity should be unrestricted, valid for anisotropic elasticity, anisotropic plasticity and mixed hardening, including the possibility of incorporating nonlinear kinematic hardening and viscoelastic effects in a simple manner for general applications. Obviously, simplicity and robustness of the algorithms are also 2

important practical requirements [23, 1]. At the time of developing a plane-stress theory and a computational algorithm for large strains, several theoretical and practical issues must be considered. First, the energetic part of the formulation should be fully hyperelastic, both in the elastic domain and in the kinematic hardening part. Therefore, formulations based on objective rates [14], tied to incrementally objective algorithms, as the Hughes-Winget [26] and Rolph-Bathe [27] algorithms, should be avoided to prevent shear stress oscillation, spurious elastic dissipation [28, 29], or the lack of weak-invariance under reference configuration changes [30]; only the logarithmic objective rate with linear elasticity (in terms of logarithmic strains) [31] fulfills Bernstein’s integrability requirements [32, 33]. In contrast, hyperelastic formulations in the context of large strain elastoplasticity [34, 35] fulfill naturally such conditions and yield simple integration algorithms. Furthermore, the elastic state variable for hyperelasticity should be obtained from Kr¨oner-Lee multiplicative decompositions (as in [34, 36, 37, 38], among others), because additive Green ansatzes and plastic metrics (as in [39, 40, 41, 42, 43, 44, 45], among others) miss a clear microstructural motivation and have been criticised for lacking weak-invariance [30] (results depend on the arbitrary reference configuration) and for loosing the elipticity properties of the elastic behavior during plastic flow [46]. Because of these reasons, we develop our formulation below in the context of multiplicative hyper-elastoplasticity. Two main approaches are traditional in the scientific community for multiplicative hyperelastoplastic formulations. The first one uses quadratic strains, it is valid for arbitrary stored energies and it may be developed quite equivalently in the intermediate or in the spatial configurations (Remark 9.2.1 in [14]); see among others [14, 36, 37, 47, 48, 49, 50]. In this approach, the computational algorithms are usually complex, specially for general anisotropic cases or for nonlinear kinematic hardening. Furthermore, the flow rule is unconventional [51]. The second approach uses logarithmic strains, see e.g. [3, 52, 53, 38], and it is frequently developed in the principal strain space [51, 54, 3]. Whereas these formulations and their associated algorithms are substantially simpler, they are restricted to elastic relations linear in logarithmic strains and to moderate elastic strains (due to the approximations employed). Both approaches use exponential mappings, either to obtain additive incremental formulations or to preserve volume. Except for the case of using principal logarithmic strains with isotropy [22, 3], to the authors’ knowledge, no plane stress projected formulation for multiplicative hyperelasto-plasticity is available and, in particular, none including nonlinear kinematic hardening. Indeed, the treatment of nonlinear kinematic hardening at large strains in the desirable hyperelastic context (see relevance of this aspect in Figs. 4-9 and 19 of [48]) is usually involved, where the kinematic state variables for such hardening are usually obtained from nested Lion multiplicative decompositions [55] of the plastic gradient into plastic energetic and plastic dissipative parts; see for example [37, 48, 56, 57, 58], among others. Then, if the 3D formulations are already complex, plane stress formulations with the constraints enforced apriori have been elusive. Recently, we have presented a novel consistent framework to deal with large strain kinematics in multiplicative elastoplasticity [59]. The proposed approach derives from a systematic use of the chain rule and introduces the concept of elastic strain rate correctors in the continuum framework. We have shown that the approach is fully consistent thermodynamically and that it can be equivalently employed using any work-conjugate pair of 3

strains and stresses; see [59] for details and for a comparison with other established formulations. Furthermore, the symmetric flow equation is conventional, the skew-symmetric flow is uncoupled, and the Mandel stress tensor plays no role. Moreover, when using logarithmic strains and generalized Kirchhoff stresses, the resulting algorithm is remarkably simple [60], even for the case of elastic and plastic anisotropy. In fact, a plain backward-Euler integration rule, without needing the exponential mapping, fully recovers the structure of the small strain algorithms and preserves the isochoric character of plastic flow. We have also developed recently a novel framework for including nonlinear kinematic hardening at large strains which is equivalent to the Ohno-Wang model (if small strains are considered) but which is fully hyperelastic in the energetic components, uses classical Kr¨oner-Lee decompositions [61, 62] (instead of Lion decompositions [55]), and does not explicitly employ the backstress concept [63]. This model and the corresponding integration algorithm are also very simple and parallel to an equivalent infinitesimal framework using classical Wilkins algorithms [64], to which simple kinematic transformations are employed. Given the success and simplicity of the novel approach, in this work we use the same concepts to develop a thermodynamically consistent theory and computational algorithm naturally constrained to the plane stress space. In the following sections, we show that, again, the resulting plane-stress constrained (rather than “projected”) formulation is fully hyperelastic (also for the energetic nonlinear hardening terms), additive, and volume preserving with a plain backward-Euler integration, keeping the simplicity of the small strain cases. Large strain kinematics are again reduced to geometric, noniterative (explicit), pre- and post-processors. The resulting formulation is remarkably parallel to both the infinitesimal one and the 3D version. We furthermore develop the consistent tangent (which also retains the same simple structure) to keep the asymptotic quadratic convergence in Newton schemes. The rest of the paper is organized as follows. First, to introduce the main concepts in the simpler infinitesimal context, we develop in Sec. 2 the plane stress-projected model for small strains using the approach to be employed at large strains, emphasizing the parallelism with the classical 3D small strains formulation found in the literature. Then, through a systematic extension of these ideas, the plane stress constrained formulation for large strains and the implicit integration algorithm are given in Secs. 3 and 4, respectively. Finally, some numerical examples are given to evaluate the performance of the proposed formulation and the congruency with the full 3D theory developed in [63]. 2. Plane-stress constrained constitutive model for small strains 2.1. Continuum formulation In this section, we formulate the plane stress theory for the infinitesimal framework. A rheological (Besseling) model as shown in Fig.1 is used, which includes N internal Prandtl (friction-spring) elements in parallel with a single Hooke (spring) element. We have shown in Ref. [63] that, without the explicit use of the backstress concept, this model is equivalent to the Ohno-Wang model, which uses multiple backstresses following each a hardening of the Armstrong-Frederick type. The tensors ε and σ are the external (i.e. measurable) infinitesimal total strains and Cauchy stresses, respectively, while the tensors εei and εpi are internal (i.e. not measurable) 4

infinitesimal elastic and plastic strains at each Prandtl device, respectively. The internal elastic strain of each Prandtl device εei can be expressed as a dependent variable of the two independent variables ε and εpi : εei (ε, εpi ) = ε − εpi

(1)

The expression εei (ε, εpi ) is general, and simply states that the internal variable εei that will be used as state variable for the stored energy, depends on the external ε and the internal εpi , the former characterizing the external evolution and the latter characterizing an internal evolution. Then, an immediate use of the chain rule brings naturally the continuum concepts of trial elastic strain rate and corrector elastic strain rate that have governed the algorithmic schemes (proposed therein as a mere algorithmic issue instead of the present continuum quantities, see e.g. [65, 14, 3] and therein references): ∂εei ∂εei e ε˙i = : ε˙ + : ε˙pi =: tr ε˙ei + ct ε˙ei (= ε˙ − ε˙pi ) (2) ∂ε ε˙p =0 ∂εpi ε=0 ˙ i

Here

tr (•)

implies trial (i.e. internal evolution frozen) and

ct (•)

implies corrector (i.e. external

(a)

(b)

Figure 1: Rheological model made of parallel Prandtl devices

evolution frozen). These rates result to be equal in the present infinitesimal case to the total strain rate and (minus) the plastic strain rate in the device, respectively. However, this identification does not longer hold in nonlinear kinematics as it will be apparent below, see a detailed analysis in [59]. From Fig. 1, the strain energy is Ψ = Ψ∞ (ε) +

N X

Ψei (εei )

= Ψ∞ (εα , ε3 (εα )) +

i=1

N X i=1

5

Ψei (εei,α , εei,3 (εei,α ))

(3)

where the sub-indices α indicate the in-plane components of the tensor and the sub-indices 3 indicate the out-of-plane components (containing the relevant non-vanishing normal component εei,33 ). Since we are considering deviatoric plasticity and the volumetric part is not affected by plastic flow, we model below all the volumetric energy as contained in the external device, i.e. in Ψ∞ . The dependencies ε3 (εα ) and εei,3 (εei,α ) are given by the plane stress condition σ3 = σ˙ 3 = 0. Thus, from Eq. (3), the strain energy rate may be obtained using the chain rule as:  N  X ∂Ψ∞ ∂Ψ∞ dε3 ∂Ψei ∂Ψei dεei,3 e e ˙ Ψ= : ε˙α + : : ε˙α + : ε˙i,α + e : e : ε˙i,α ∂ε ∂ε dεα ∂εe ∂εi,3 dεi,α | α {z 3 } i=1 | i,α {z } ˙ ∞ =(dΨ∞ /dεα ):ε˙α Ψ

(4)

˙e ˙ e +ct Ψ ˙ e =tr Ψ Ψ i i i

The dissipation rate DP constrained to the plane stress subspace can be determined from ˙ according to the Clausiusthe external stress power P and the total strain energy rate Ψ Duhem inequality as ˙ ≥0 DP = P − Ψ (5) Due to the plane stress conditions (σ3 = 0), since there are no out-of-plane loads, the external stress power can be determined as P = σ : ε˙ = σα : ε˙α , the latter scalar product containing only in-plane (•)α components. Two cases can be analyzed for the inequality in Eq. (5). In the first case, the rate of the internal variable vanishes even if ε˙ 6= 0, i.e. ε˙pi = 0 for all i = 1, ..., N , so no dissipation takes place in any Prandtl device. Then, from Eq. (1) ε˙ei = tr ε˙ei = ε˙ = {ε˙α , ε˙3 }, ∀i and from Eq. (5) we have  N X ˙ e ˙ = σα : ε˙α − Ψ ˙ kin + Ψ D = σα : ε˙α − Ψ i P

tr

i=1

ε˙pi =0

 =0

if all

ε˙pi = 0

(6)

tr Ψ ˙

˙ ˙p where ≡ Ψ| εi =0,∀i is the “trial” rate (plastic, internal evolution frozen) and which for the case ε˙α 6= 0 requires dΨ dε σ α : ε˙α = : : ε˙α (7) dε ε˙p =0,∀i dεα i

For general ε˙α 6= 0 this means that N dΨ dε dΨ dΨ∞ X dΨei ∂εei,α : = = + : σα = dε ε˙p =0,∀i dεα dεα ε˙p =0,∀i dεα dεei,α ∂εα ε˙p =0 i=1 i i i ! N e e ∂εi,α ∂Ψ∞ ∂Ψ∞ dε3 X ∂Ψei ∂Ψei dεi,3 = + : + + : : ∂εα ∂ε3 dεα ∂εei,α ∂εei,3 dεei,α ∂εα ε˙p =0 i | {z } i=1 | {z } =: σ∞,α |e =: σi,α N X ∂εei,α |e = σ∞,α + σi,α : ∂εα ε˙p =0 i=1

i

6

(8)

|e

Note that each in-plane stress tensor (σi,α or σ∞,α ) is composed of two addends. The first one is the partial derivative of the corresponding strain energy respect to the in-plane stresses (i.e. the component typical in 3D analysis). The second addend is the correction due to the plane stress constraint; i.e. an in-plane strain produces an out-of-plane strain which also contributes to the stored energy (but not to the external power because out-of-plane stresses vanish). These expressions reveal the fact that we are considering general 3D strain energy functions. In-plane, already constrained energy expressions could also be used and some savings could be obtained, but they are usually not so readily available in general large strain analyses. If in-plane strains remain constant (ε˙α = 0), for the constrained problem the rate of the stored energy cannot change for changes in the out-of-plane strains ε˙3 6= 0, ˙ i.e. (∂ Ψ/∂ε 3 ) : ε˙ 3 = 0 (or alternatively the rate of the out-of-plane ∂Ψ∞ /∂ε3 stress should vanish), see Eq. (6). Then   ˙∞ ∂Ψ d ∂Ψ∞ ∂ 2 Ψ∞ ∂ 2 Ψ∞ : ε˙3 ≡ : ε˙3 = 0 = ε˙3 : : ε˙α + ε˙3 : : ε˙3 (9) ∂ε3 dt ∂ε3 ∂ε3 ∂εα ∂ε3 ∂ε3 and also by Eq. (1), when ε˙pi = 0, we have tr ε˙ei = ε˙ and # " 2 Ψe 2 Ψe ˙ e ∂εei,3 ∂εei,3 ∂Ψ ∂ ∂ e e i i i ˙ ˙ : : ε + : ε : ε ˙ = 0 = : : ε˙3 3 i,α i,3 ∂εei,3 ∂ε3 ∂εei,3 ∂εei,α ∂εei,3 ∂εei,3 ∂ε3

(10)

Then, we have that since strain rates ε˙3 may take any general value, the dependencies dε3 /dεα and dεei,3 /dεei,α must take the relations  2 −1 ∂ Ψ∞ ∂ 2 Ψ∞ dε3 =− : (11) dεα ∂ε3 ∂ε3 ∂ε3 ∂εα !−1 dεei,3 ∂ 2 Ψei ∂ 2 Ψei = − : (12) dεei,α ∂εei,3 ∂εei,3 ∂εei,3 ∂εei,α The in-plane internal stress rate for each Prandtl device is: tr |e σ˙ i,α

d = dt  =

∂Ψei ∂Ψei dεei,3 + : ∂εei,α ∂εei,3 dεei,α

∂ 2 Ψei ∂ 2 Ψei − ∂εei,α ∂εei,α ∂εei,α ∂εei,3

|

!

dεei,3 ∂ 2 Ψei ∂ 2 Ψei + : = ∂εei,α ∂εei,α ∂εei,α ∂εei,3 dεei,α ε˙pi =0  !−1 2 e 2 e ∂ Ψi ∂ Ψi  tr e : : : ε˙i,α e e ∂εi,3 ∂εi,3 ∂εei,3 ∂εei,α {z } |e ¯ Ci,α

|e

¯ : tr ε˙e =C i,α i,α

! : tr ε˙ei,α

(13)

where performing the time derivative in parenthesis we have used the identity in Eq. (10) to ¯ |e is defined as the internal plane stress-projected elastoplastic cancel some terms and where C i,α tangent moduli. The trial in-plane internal stresses can be determined as tr

|e

|e

σi,α = t σi,α +

τ ∈[t,tr]

7

¯ |e : (tr εe − t εe ) C i,α i,α i,α

(14)

where τ represents the “time” at which the tangent is evaluated (typically τ = tr for backward-Euler). Obviously, σi,3 = 0. For the usual cases of metal plasticity, the tangent moduli are constant (e.g. from the plane stress linear elasticity matrix). If the stored energy ¯ ∞ (εα ), then we immediately have is already written in a plane-stress constrained form, e.g. Ψ ¯ for example σ∞,α = dΨ∞ (εα )/dεα , and so on. In order to be illustrative, as a simple (but usual) example, consider the stored energy as a function of the deviatoric principal strain components εed i,k , k = 1, ..., 3, namely ed ed 2 ed 2 ed 2 Ψei (εei,1d , εed i,2 , εi,3 ) := µ(εi,1 ) + µ(εi,2 ) + µ(εi,3 )

In the main 3 × 3 box of  2  " ed T d Ψ dεi = dεei dεei dεei  2 − 31 3 2 = − 13 3 1 − 3 − 31

principal directions, we have # dεed d2 Ψ i : ed ed : dεei dεi dεi   2 − 13 2µ − 13 3 2  − 1 − 31   2µ 3 3 2 1 2µ − 3 − 31 3

(15)

  4  − 13 − 23 µ − 32 µ 3µ − 13  = − 32 µ 43 µ − 32 µ (16) 2 − 32 µ − 23 µ 34 µ 3

and from Eq. (13), taking into account that in these principal directions representation the third dimension is the out-of-plane one, we have   h i  4 µ − 2 µ − 2 µ     µ −µ −1 |e 4 2 2 3 3 3 ¯ µ −3µ −3µ = (17) − Ci,α = − 23 µ 43 µ − 23 µ 3 −µ µ |e ¯ |e : tr εe and tr σi,3 = 0; note that also tr(tr σ |e ) = tr(tr σ |e ) = 0, so in this case tr σi,α = C i,α i,α i i,α which means that both the complete and the in-plane stresses are deviatoric. Of course   tr e  1 tr εe = 1 ε , so the strains do not need to be isochoric. Remarkably, using isochoric i,α i,3 2 2 strain energy functions of the type Eq. (15) in the internal Prandtl devices, we obtain |e |e |ed σi = {σi,α ≡ σi,α , σi,3 = 0} that are both deviatoric and constrained to the plane stress space, regardless of εei having any volumetric component. This is central to the simplicity of the resulting formulation. Consider now the case in which the system is mechanically isolated, ε˙ = 0, and only an internal evolution of the independent variable takes place ε˙pi 6= 0. The dissipation inequality of Eq. (5) must be positive due to the burgeoning of the non-zero plastic deformation ! N N N e e e X X X dε ∂Ψ ∂Ψ ct i,3 |e P e ct e i i ˙ ˙ D = − Ψ =− Ψi = − + e : e : ε˙i,α =: − σi,α : ct ε˙ei,α e ˙ ∂εi,α ∂εi,3 dεi,α ε=0 i=1 i=1 i=1 (18) In our rheological model, we have defined an equivalent stress-like variable ki and a strainlike rate γ˙ i such that the dissipated power of each Prandtl device is ki γ˙ i . This definition is ˙ i (γi ) =: ki γ˙ i is equivalent to state that the dissipation potential in each device is Φi (γi ) so Φ the dissipation rate. Therefore for each device we have the following relation: ct Φ˙ i (γi ) = − Ψ˙ ei

8

(19)

or

|e δ˙i := −σi,α : ct ε˙ei,α − ki γ˙ i = 0

(20)

where δ˙i is the balance in the dissipation, i.e. the difference between stored energy loss rate and dissipation rate. In order to determine ct ε˙ei,α , without loss of generality, let us assume the following form for ct ε˙e,α : r 3 ct ˆ i,α ε˙e,α = − γ˙ i n (21) 2 p ˆ i,α is a yet undetermined in-plane unitary tensor and the factor 3/2 is introduced in which n for convenience in using typical (familiar to the reader) uniaxial quantities. Using Eq. (21), each addend in Eq. (18) is ! r r 3 3 |e |e ˆ i,α γ˙ i = ki γ˙ i =⇒ ˆ i,α − ki γ˙ i = 0 (22) DiP = σ :n σ :n 2 i,α 2 i,α | {z } q =: 32 fi so either γ˙ i = 0 or fi = 0. The principle of maximum dissipation implies that the in-plane ˆ i,α is such that the dissipation is maximum, i.e. taking the normalizing constraint direction n ˆ i,α : n ˆ i,α − 1 = 0 C=n d(DP − λC) =0⇒ ˆ i,α dn

r

3 |e ˆ i,α = 0 ⇒ n ˆ i,α = σ − 2λn 2 i,α

r

|e

σi,α 2 1 |e σi,α = |e 3 2λ ||σi,α ||

(23) |e

where the value of λ has been obtained from the constraint C = 0; note that since σi,3 = 0, |e

|e

we have that ||σi,α || ≡ ||σi ||. Then |e

ˆ i,α − fi = σi,α : n

q

2 3 ki

=0

(24)

is the yield function for the plane-stress problem, which must hold when γ˙ i 6= 0. From Eqs. (21) and (23), the flow rule is then given as r ct e ε˙i,α

=−

|e

σi,α 3 γ˙ i |e 2 ||σ ||

(25)

i,α

Note that in Eq. (24) the stress only contains the in-plane components with the plane-stress ˆ i,α is also defined in that plane with the same constraint already enforced, and the normal n constraint. Hence, the reader should not be mislead by the (remarkable, advantageous and pursued) same appearance of Eqs. (24) and (25) to those of the 3D theory. They are not simply obtained neglecting the out-of-plane components or “projecting” them (even though we keep the nomenclature of “projected” formulation). Indeed, these are in contrast to those found, for example, in Eq. (2.4.19) of Ref. [14], which are a phenomenological projection of those given in Eq. (2.3.4) of Ref. [14]. Noteworthy, Eqs. (24) and (25) have been derived systematically from the plane stress requirements in the general hyperelastic framework which 9

will be used in the large strain context below. For instance, in Eq. (2.4.10) of Ref. [14], both ˜ is neither the stress σ and the overstress η are non-deviatoric, and hence, the backstress β deviatoric. Then, η needs to be projected to the plane stress space to be used both in the flow |e rule and in the yield stress. In contrast, in our Eqs. (24) and (25), σi,α is purely deviatoric and in-plane by construction, conditions derived from the model structure, thermodynamic consistency, and from the principle of maximum dissipation. Of course, the similarity of Eqs. (24) and (25) with the classical 3D equations is a very interesting asset of the present formulation. By applying the chain rule of differentiation, the rate of the yield function may be obtained as f˙i = f˙i |ε˙pi =0 + f˙i |ε=0 = tr f˙i + ct f˙i (26) ˙ For the first addend, which stands for the trial elastic evolution, we have f˙i |ε˙pi =0 =

tr

|e ¯ |e : tr ε˙e ˆ i,α : tr σ˙ i,α = n ˆ i,α : C f˙i = n i,α i,α

(27)

For the second addend in Eq. (26), which stands for a corrector elastic evolution, we have q q |e ct ˙ ct |e ct e 2 ct ˙ 2 ct ˙ ¯ ˙ ˆ ˙ ˆ σ − f˙i |ε=0 = f = n : ε − k = n : C : ˙ i i,α i i,α i,α i,α i,α 3 3 ki q q 2 0 ¯ |e : n ˆ i,α : C ˙i (28) = − 32 γ˙ i n i,α ˆ i,α − 3 ki γ where we have included the possibility of isotropic hardening by a function ki (γi ). From the consistency condition, γ˙ i > 0 requires f˙i = 0 at any time t, and therefore we can factor-out γ˙ i and obtain the following familiar relation r γ˙ i =

¯ |e : tr ε˙e ˆ i,α : C 2 n i,α i,α |e 3n ¯ :n ˆ i,α : C ˆ i,α + 2 k 0 i,α

if ε˙pi 6= 0

(29)

3 i

Combining Eq. (29) with the associative flow rule Eq. (25), we can establish the consistent coupling between the internal corrector elastic strain rate ct ε˙ei,α and the internal trial elastic strain rate tr ε˙ei,α ! ! ¯ |e ) ¯ |e : tr ε˙e ˆ i,α ⊗ (n ˆ i,α : C ˆ i,α : C n n i,α i,α i,α ct e ˆ i,α = − ε˙i,α = − n : tr ε˙ei,α (30) |e |e 2 0 2 0 ¯ ¯ ˆ i,α : Ci,α : n ˆ i,α + 3 ki ˆ i,α : Ci,α : n ˆ i,α + 3 ki n n If we now substitute Eq. (30) into Eq. (2), we have the relation between the internal elastic strain rate and the trial elastic strain rate during plastic flow ! ¯ |e ) ˆ ˆ n ⊗ ( n : C i,α i,α i,α : tr ε˙ei,α if γ˙ i > 0 (31) ε˙ei,α = ¯IS − |e 2 0 ¯ ˆ ˆ ni,α : Ci,α : ni,α + 3 ki   being ¯IS αβγ = 12 (δαγ δβ + δα δβγ ) the in-plane four-order symmetric identity tensor, with α, β, γ,  = 1, 2, being the in-plane indices. Note that the appearance of all these expressions is the same as in the 3D theory but considering plane-constrained entities. Hence, further details are omitted. 10

2.2. Integration algorithm The derivation of the integration algorithm, based on the closest point projection method [66], is also straightforward, and in some sense, it is similar to the additive algorithms in strain space [67]. The standard backward Euler method is used in the plastic corrector step for each device i to compute the evolution of the plastic strain, so the two residuals (one related to the elastic strain tensor and the other one related to the yield condition) can be set as

 Ri (Ei ) :=

ρi,α fi



  q   t+∆t εe − tr εe + 3 ∆γi t+∆t n  ˆ i,α 0 i,α i,α 2 q = → 2 t+∆t 0  t+∆t σ |e : t+∆t n  ˆ i,α − ki i,α

(32)

3

with Ei being the vector of variables to be determined for device i, defined as Ei := [εei,α , γi ]T . The solution of the above system can be efficiently performed by the Newton-Raphson algorithm; in fact the highly efficient Wilkins algorithm [64] may be used. If needed, the out-of-plane internal strain εei,3 may be updated at each iteration using Eq. (12). Since the steps to derive the iterative algorithm are standard, we omit further details. 3. Plane-stress constrained constitutive model: large strains formulation 3.1. Kinematics and logarithmic strain rate tensors In finite strain plasticity, motivated from crystal plasticity, the Kr¨oner-Lee multiplicative decomposition of the deformation gradient F into an elastic part F e and a plastic part F p is widely accepted. We use here the ideas of the foregoing rheological model presented for the small strains formulation, so a multiplicative split at each Prandtl device i is employed F = Fie Fip

(33)

The total Green-Lagrange strain tensor, the plastic Green-Lagrange strain tensor in the reference configuration, and the elastic Green-Lagrange strain tensor in the intermediate configuration are defined as T

T

A := 12 (F T F − I), Api := 21 (Fip Fip − I), and Aei := 12 (Fie Fie − I)

(34)

where I is the second order identity tensor, with components being the Kronecker deltas, [I]ij = δij . Following the ideas introduced in the small strains formulation, we consider the dependent internal elastic variable Aei to be a function of two independent variables: the external strain tensor A and the internal plastic deformation gradient tensor Fip , as follows—the reader may verify that this identity holds using the definitions in Eq. (34) −T

Aei (A, Api ) = Fip

(A − Api )Fip

−1

= Fip

−T

Fip

−T

: (A − Api )

(35)

where the operator is the mixed dyadic product between two second order tensors, i.e. in components (Y Z)ijkl = Yik Zjl . Note that Eq. (35) simply states that for being additive, 11

Box 1 Implicit integration algorithm for small strains Require: Given the state variables at time t and the total strain at time t + ∆t. 1: Elastic predictor: evaluate the trial state tr εe = t+∆t ε − t εp tr εe = t εe α i,α i,3 i,3 i,α and tr σ |e i,α

|e

|e

= t σi,α + tr Ci,α : (tr εei,α − t εei,α ); 2: Check yield condition and test for plastic loading |e tr ki := t ki ; tr fi = fi (tr σi,α , tr ki ) IF tr fi ≤ 0 THEN Elastic step: set t+∆t (•) = tr (•) and EXIT ELSE Plastic step: Proceed to step 3 ENDIF 3: Plastic return mapping (i) Plastic corrector: Solve the system of equation for {t+∆t εei,α , t+∆t γi }   q t+∆t εe − tr εe + 3 ∆γi t+∆t n    ˆ i,α 0 i,α i,α 2 q = |e t+∆t t+∆t 2 0 t+∆t  ˆ i,α − ki  σ : n i,α

3 t+∆t εe,(0) i,α

(0)

with initial values = tr εei,α , t+∆t γi = tr γi and (ii) Update the internal stresses and plastic tangent moduli t+∆t

|e

|e

σi,α = t σi,α +

t+∆t ¯

C=

t+∆t ¯ |e Ci,α

t+∆t ¯

C∞ +

4:

EXIT

12

: (t+∆t εei,α − t εei,α )

N X t+∆t i=1

t+∆t (0) ki

¯ |e : ¯IS C i,α

= ki (t γi )

strains need to be placed in the same configuration (A and Api lie in the reference configuration, whereas Aei lies in the intermediate configuration of device i). Whereas Green-Lagrange strains are accessible immediately from the deformation gradients, they are not advantageous for developing simple formulations. In contrast, the natural properties of logarithmic strains bring remarkably simple theoretical and computational frameworks in elasto-plasticity; hence they are our selection. Since one-to-one relations Eie = Eie (Aei ) and E = E(A) must hold (deformation states T are unique), where Eie = 21 ln(Fie Fie ) and E = 12 ln(F T F ) are the elastic and total “material” (or reference) logarithmic strain tensors in their respective configurations, we can also formulate the general dependence Eie (E, Fip ) which generalizes Eq. (1) and Eq.(35) above to the logarithmic strain space. Hence, in a parallel form as in Eq. (2), by immediate application of the chain rule, we can decompose the elastic logarithmic strain rate tensor into the addition of two partial contributions: ∂Eie ∂Eie e ˙ ˙ :E+ : F˙ip = E˙ ie |F˙ p =0 + E˙ ie |E=0 Ei = (36) ˙ i ∂E F˙ p =0 ∂Fip E=0 ˙ i

Using the terminology employed typically in the operator split of computational inelasticity (there defined in the algorithmic framework but here used in the continuum setting), we can write E˙ e = E˙ e | ˙ p + E˙ e | ˙ = tr E˙ e + ct E˙ e (37) i

i Fi =0

i E=0

i

i

where tr E˙ ie is the trial contribution to the rate of the total elastic logarithmic strain and E˙ ie depends on the total logarithmic strain rate E˙ only (i.e. the plastic evolution given by Fip is frozen), and ct E˙ ie is the plastic corrector contribution to the total elastic logarithmic strain rate E˙ ie and depends on the total plastic deformation gradient rate F˙ip only (i.e. the external deformation given by E is frozen). Note that in contrast to what happened in Eq. (2), here E˙ ie 6= E˙ − E˙ ip (if E p = 21 ln(F pT F p ) or similarly defined) except in the case of proportional loads. Indeed, Eq. (37) is an immediate consequence of the chain rule, but ansatzes of the type E˙ ie = E˙ − E˙ ip are incompatible with the multiplicative decomposition and present several fundamental issues (loss of elipticity, lack of weak invariance), as explained in the Introduction Section. Furthermore, we remark that we will never define nor use a measure of the type Eip in our formulation; only the “elastic” strains Eie and the dissipated energy rate ki γ˙i are relevant. 3.2. Dissipation inequality and flow rule in terms of corrector elastic strain rates By using work-conjugacy, the stress power per reference volume may be written in terms of the (“material”) generalized Kirchhoff stress tensor T and its work-conjugate logarithmic strains E (see Ref. [68]) as P = τ : d = S : A˙ = T : E˙ (38) where τ are the Kirchhoff stresses, d is the deformation rate tensor, and S are the 2nd Piola-Kirchhoff stresses. For the plane stress state, due to the fact that the rotation tensor is in-plane (so the third direction remains invariant), and that there are no loads in that axis, we can write in a similar way as in the infinitesimal case P = T : E˙ = Tα : E˙ α . 13

Following also the same framework, the rheological model in Fig. 1 motivates the additive e decomposition of the stored energy Ψ into the internal elastic part Ψ Pi Nof Ne Prandtl devices and a part associated with the limiting hardening Ψ∞ , i.e. Ψ = Ψ∞ + i=1 Ψi . The kinematic stored energy function Ψ∞ (E) is a function of the total logarithmic strain E, and includes the volumetric part, while the internal elastic stored energy Ψei (Eie ) function is a function of the internal elastic logarithmic strain Eie . Then —cf. Eq. (3) Ψ = Ψ∞ (E) +

N X

Ψei (Eie ) = Ψ∞ (Eα , E3 (Eα )) +

i=1

N X

e e e Ψei (Ei,α , Ei,3 (Ei,α ))

(39)

i=1

e ˙ =Ψ ˙ ∞ (Eα , E3 ) + PN Ψ ˙e e and its rate form is Ψ i=1 i (Ei,α , Ei,3 ). As with the infinitesimal strain formulation, the stored energy rate of each device i may be expressed in terms of the trial and corrector contributions

˙ e (E e , E e ) = tr Ψ ˙ e (E e , E e ) + ct Ψ ˙ e (E e , E e ) Ψ i i,α i,3 i i,α i,3 i i,α i,3

(40)

According to the Clausius-Duhem inequality, the dissipation for the plane stress problem can be determined as ˙ DP = P − Ψ " ˙ ∞ (Eα , E3 ) + = Tα : E˙ α − Ψ

N X

# tr

(

e e ˙ ei (Ei,α Ψ , Ei,3 )

+

ct

e e ˙ ei (Ei,α Ψ , Ei,3 ))

≥0

(41)

i=1

Consider the first case for which F˙ip = 0 for all Prandtl devices, i.e. there is no dissipation, ˙ Then, ct Ψ ˙ e = 0 and the dissipation and therefore we have E˙ ie ≡ tr E˙ ie = ∂Eie /∂E|F˙ p =0 : E. i i inequality is ˙ ∞ (Eα , E3 ) − D = Tα : E˙ α − Ψ P

N X

tr

˙ e (E e , E e ) = 0 Ψ i i,α i,3

if F˙ip = 0 for all i

(42)

i=1

The external stored energy rate is determined as   ∂Ψ∞ ∂Ψ∞ dE3 dΨ∞ ˙ ˙ : Eα = + : : E˙ α := T∞,α : E˙ α Ψ∞ (Eα , E3 ) = dEα ∂Eα ∂E3 dEα {z } | T∞,α and the ith trial internal stored energy rate can also be determined as e ∂Ei,α dΨei tr ˙ e e e : E˙ α Ψi (Ei,α , Ei,3 ) = e : ∂E p dEi,α α F˙ =0 i ! e e e e dEi,3 ∂Ei,α ∂Ψi ∂Ψi ˙ = : e + ∂E e : dE e ˙ p : Eα ∂Ei,α ∂E α i,3 i,α Fi =0 | {z } |e := T i,α e ∂Ei,α |e = Ti,α : : E˙ α ∂Eα F˙ p =0 i

14

(43)

(44) (45)

(46)

so for general cases of non-dissipative deformation, Eq. (42) requires N e X ∂Ei,α = T∞,α + Ti,α : ∂Eα F˙ p =0

(47)

e e , E e (E e )) dEi,3 dΨei (Ei,α ∂Ψei ∂Ψei i,α i,α ≡ + : e e e e dEi,α ∂Ei,α ∂Ei,3 dEi,α

(48)

Tα = T∞,α +

N X

|e Ti,α

i=1

i=1

i

in which we defined |e

Ti,α := and

Ti,α :=

|e Ti,α

e ∂Ei,α : ∂Eα F˙ p =0

(49)

i

are the ith internal independent stresses in the reference configuration. Aside, following the same steps as in the infinitesimal formulation, we get  2 −1 dE3 ∂ Ψ∞ ∂ 2 Ψ∞ =− : (50) dEα ∂E3 ∂E3 ∂E3 ∂Eα !−1 e dEi,3 ∂ 2 Ψei ∂ 2 Ψei = − : (51) e e ∂E e e ∂E e dEi,α ∂Ei,3 ∂Ei,3 i,3 i,α so, for example, following again the same steps to arrive at Eq. (13), we have e , Ee ) d2 Ψei (Ei,α i,3 e ¯ |e : E˙ e : E˙ i,α =: A i,α i,α e dE e dEi,α i,α   !−1 2 e 2 e 2 e 2 e ∂ Ψi ∂ Ψi ∂ Ψi ∂ Ψi e  : E˙ i,α = − : : e e e e ∂E e e e e ∂Ei,α ∂Ei,α ∂Ei,3 ∂Ei,α ∂Ei,3 ∂Ei,3 ∂Ei,3 i,α

|e T˙i,α =

(52)

(53)

Furthermore, the same stored energy of Eq. (15) and reduced constitutive tensor in Eq. (17) may be used if written in terms of logarithmic strains (note that infinitesimal strains and logarithmic strains share the same structure, e.g. deviatoric and volumetric operators). This type of constitutive equation, and its nonlinear and anisotropic (linear and nonlinear) extensions, have already been used in polymers [69], foams and auxetic materials [70], and in anisotropic soft tissues [71, 72], including finite nonlinear viscoelasticity [73]. In particular, if we assume the classical Valanis-Landel decomposition employed in isotropic polymers, which is the nonlinear generalization of Eq. (15), we have ed ed ed ) + ωi (Ei,2 ) + ωi (Ei,3 ) Ψei (Eie ) = ωi (Ei,1

(54)

ed is the principal component k of E ed , and ω (•) where Eied is the deviatoric part of Eie , Ei,k i i is a function. Then, the main 3 × 3 box is #  "  h i 2 Ψe d µ ¯i −¯ µi |e i ¯ = (55) Ai,α := e dE e −¯ µi µ ¯i dEi,α i,α

15

where µ ¯i =

ed )ω 00 (E ed ) + ω 00 (E ed )ω 00 (E ed ) + ω 00 (E ed )ω 00 (E ed ) ωi00 (Ei,1 i i,2 i i,1 i i,3 i i,2 i i,3 ed ) + ω 00 (E ed ) + 4ω 00 (E ed ) ωi00 (Ei,1 i i,2 i i,3

(56)

ed = −E ed − E ed , and we used the notation and Ei,3 i,1 i,2 ed ωi00 (Ei,k )

d2 ωi (E d ) = d(E d )2 E d =E ed

(57)

i,k

|e which for the linear case results in µ ¯i = µi ≡ ωi00 /2. It is obvious that T˙i,α is deviatoric and |e |e e , and recall that we imposed T˙i,3 = 0, so Ti is always in-plane and deviatoric. Of course Ei,3 e its relevant component Ei,3 , may be obtained from Eq. (51) in rate form as e E˙ i,3 =−

ed ) − ω 00 (E ed )]E ˙ e + [2ω 00 (E ed ) − ω 00 (E ed )]E˙ e + 2ω 00 (E ed )(E˙ e + E˙ e ) [2ωi00 (Ei,1 i i,2 i,1 i i,2 i i,1 i,2 i i,3 i,1 i,2 ed ) + ω 00 (E ed ) + 4ω 00 (E ed ) ωi00 (Ei,1 i i,2 i i,3

(58) For the second case, E˙ = 0, so the external deformation given by E is frozen and there may be a plastic evolution for each device F˙ip 6= 0. In this case, E˙ ie = ct E˙ ie = ∂Eie /∂Fip |E=0 : F˙ip , ˙ and the dissipation inequality, Eq. (41), simplifies to ! N N e X X dEi,3 ∂Ψei ∂Ψei |e P ct ˙ e e D =− Ψi = − : ct E˙ ie = −Ti,α : ct E˙ i,α > 0 for F˙ip 6= 0 e + ∂E e : dE e ∂Ei,α i,3 i,α i=1 i=1 (59) In contrast to other theories, it is expressed in terms of symmetric tensors of elastic nature lying in the intermediate configuration, namely the corrector contribution to the elastic loge and its power-conjugate generalized Kirchhoff stress tensor T |e . arithmic strain rate E˙ i,α i,α Note that the Mandel stress tensor plays no role and the formulation is uncoupled from any skew-symmetric contribution; see analysis in [59]. 3.3. Yield function The dissipation inequality obtained in Eq. (59) should be positive for all possible motions according to the second law of the thermodynamics, which imposes the restrictions on the possible forms of the plastic evolution. Herein, following the same arguments that arrived to Eq. (20), we arrive at the following flow rule q ct e ˆ i,α E˙ i,α = − 32 γ˙ i N (60) ˆ i,α = T |e /||T |e || is the in-plane and deviatoric associative flow direction. Then, where N i,α i,α if we now define ki as a Kirchhoff-stress-like uniaxial yield stress for the device and γ˙ i as the work-conjugate uniaxial plastic strain rate, following once more the same steps as in the infinitesimal case, we have DP =

N q X

|e 3 ˙ i Ti 2γ

ˆ i,α ≡ :N

N X i=1

i=1

16

ki γ˙ i > 0

if γ˙ i > 0

(61)

which results in q |e ˆ i,α − ki )γ˙ i = 0 ⇐⇒ f¯i γ˙ i = 0 ( 32 Ti : N | {z } q

if γ˙ i > 0

(62)

3¯ 2 fi

and the yield function f¯i is identified as a function of the in-plane part (•)α of the internal generalized Kirchhoff stress tensor for the plane stress problem, that is q |e |e ˆ i,α − 2 ki = 0 if γ˙ i > 0 f¯i (Ti,α , ki ) = Ti,α : N (63) 3 To derive the algorithm below, the substitution of q ˆ i,α = 3 T |e /ki N 2 i,α

(64)

in both the flow rule and the yield function gives simpler manipulations (specially when no isotropic hardening takes place), i.e. in the derivations we will use the also classical form— note that fi 6= f¯i but fi = 0 ⇐⇒ f¯i = 0 |e

|e

|e

fi (Ti,α , ki ) = Ti,α : Ti,α − 23 ki2 = 0 and ct

if γ˙ i > 0

γ˙ i |e |e E˙ i,α = − 32 Ti,α ki

(65)

(66)

3.4. Linearization of the continuum theory: continuum tangent moduli in the plane-stress subspace From Eq. (47), the in-plane total stresses Tα include two contributions: one from the external kinematic stresses T∞,α and one from the N internal isotropic stresses Ti,α . In rate form, we have the following relation T˙α = T˙∞,α +

N X

T˙i,α

(67)

i=1

The stresses T∞,α are determined from the strain energy function Ψ∞ (Eα , E3 ) once the logarithmic strains E are obtained from the deformation gradient tensor F . Hence, the rate form of T∞,α in the reference configuration is determined as   2 ∂ 2 Ψ∞ dE3 ∂ Ψ∞ ¯ ∞,α : E˙ α ˙ T∞,α = + : : E˙ α =: A (68) ∂Eα ∂Eα ∂Eα ∂E3 dEα ¯ ∞,α contains the external kinematic tangent moduli for the plane-stress problem. By where A using the dependence in Eq. (50), these moduli are ¯ ∞,α := A

∂ 2 Ψ∞ ∂ 2 Ψ∞ − : ∂Eα ∂Eα ∂Eα ∂E3 17



∂ 2 Ψ∞ ∂E3 ∂E3

−1 :

∂ 2 Ψkin ∂E3 ∂Eα

(69)

As mentioned previously, the stresses tensor Ti,α represents the pull-back operation ap|e plied to Ti,α from the intermediate configuration to the reference configuration. Therefore, in order to determine the rate form of Ti,α in the reference configuration, we have to determine |e the rate form of Ti,α in the intermediate configuration of device i, i.e. " # e 2 Ψe 2 Ψe dE ∂ ∂ i,3 |e e i i ¯ |e : E˙ e T˙i,α = : E˙ i,α =: A (70) i,α i,α e ∂E e + ∂E e ∂E e : dE e ∂Ei,α i,α i,α i,α i,3 |e

¯ where we defined A i,α to be the in-plane hyperelastic constitutive tensor for logarithmic strains in the intermediate configuration of device i. By conservation principles, Eq. (62) must hold at all times, and since fi = 0 ⇐⇒ f¯i = 0, we require d(fi γ˙ i ) = 0 ⇐⇒ f˙i γ˙ i + fi γ¨i = 0 =⇒ f˙i = 0 if γ˙ i > 0 (71) dt The rate of the yield function may be formulated in two parts –as in Eq. (26)– as follows f˙i = f˙i |F˙ p =0 + f˙i |E=0 = ˙ i

where 1 ˙ p 2 fi |F˙i =0

=

1 tr ˙ 2 fi

|e

= Ti,α :

tr

tr

f˙i +

ct

f˙i

(72)

|e |e ¯ |e : T˙i,α = Ti,α : A i,α

tr

e E˙ i,α

(73)

and 1 ˙ ˙ 2 fi |E=0

=

1 ct ˙ 2 fi

|e |e |e ¯ |e : = Ti,α : ct T˙i,α − 23 ki k˙ i γ˙ = Ti,α : A i,α 1 |e ¯ |e : T |e − 4 k 2 k 0 ) = γ˙ 2 (−Ti,α : A i,α i,α 9 i i k 3 i

ct

E˙ e,α − 23 ki k˙ i γ˙ (74)

Then, from Condition (71), we have the following relation for the plastic consistency parameter ! |e ¯ |e : tr E˙ e T : A i,α i,α i,α γ˙ i = 32 ki if F˙ip 6= 0 (75) |e |e |e 4 2 0 ¯ Ti,α : Ai,α : Ti,α + 9 ki ki Substituting Eq. (75) into Eq. (66), we have ˙e |˙ E i,α E=0

=

ct

e E˙ i,α

=−

|e ¯ |e : T |e ) Ti,α ⊗ (A i,α i,α |e Ti,α

¯ |e : T |e + 4 k 2 k 0 :A i,α i,α 9 i i

! :

Now, using Eq. (76) into Eq. (37), we arrive at " # |e ¯ |e : T |e ) Ti,α ⊗ (A i,α i,α e S ¯ ˙ E : i,α = I − |e ¯ |e : T |e + 4 k 2 k 0 Ti,α : A i,α i,α 9 i i

tr

e E˙ i,α

tr

e E˙ i,α

if γ˙ i > 0

(76)

(77)

or alternatively, reverting the change in Eq. (64) we find an expression with the same aspect as in Eq. (30) " # ¯ |e : N ˆ i,α ⊗ (A ˆ i,α ) N i,α e ˙ e = ¯IS − E : tr E˙ i,α (78) i,α |e 2 0 ¯ ˆ ˆ Ni,α : Ai,α : Ni,α + 3 ki 18

It can be observed that the term in brackets [•] above maps the trial rate into the final one. Hence, the in-plane rate of internal stress in Eq. (70) can be rewritten as " # ¯ |e : T |e ) ⊗ (A ¯ |e : T |e ) (A i,α i,α i,α i,α |e |e e ¯ − ¯ |ep : tr E˙ e T˙i,α = A : tr E˙ i,α =A (79) i,α i,α i,α |e |e |e 4 2 0 ¯ Ti,α : Ai,α : Ti,α + 9 ki ki ¯ |ep as the plane-stress constrained version of the internal continuum where we interpret A i,α elastoplastic tangent tensor, lying in the actual intermediate configuration. Using the relation |e between Ti,α and Ti,α established in Eq. (49), the rate of the internal stresses tensor Ti,α in the reference configuration can be obtained by applying the chain rule of differentiation as ! e e ∂E ∂E d i,α i,α |e |e (80) + Ti,α : T˙i,α = T˙i,α : ∂Eα F˙ p =0 dt ∂Eα F˙ p =0 i

i

We note that the second term in Eq. (80) is the typical convective geometric term needed to obtain the relation of T˙i,α with E˙α , which is difficult to obtain due to the continuous change of the intermediate configuration. However, during the algorithmic stress integration the trial configuration will remain fixed as in the Lie derivative, so only the first addend is of interest and the following relation can be established for the plane stress subspace e e T e ∂Ei,α ∂Ei,α ∂Ei,α |e |ep ¯ i,α : E˙ α (81) ¯ ˙ ˙ Ti,α |F˙ p =0 = Ti,α : = : E˙ α = A : Ai,α : i ∂Eα F˙ p =0 ∂Eα F˙ p =0 ∂Eα F˙ p =0 i

i

i

¯ i,α as the plane-stress version of the continuum elastoplastic moduli in the where we define A |ep reference configuration by means of the pull-back operation over its internal counterpart Ai,α at the intermediate configuration e T e ∂Ei,α ∂Ei,α |ep ¯ ¯ Ai,α = : Ai,α : (82) ∂Eα F˙ p =0 ∂Eα F˙ p =0 i

i

The complete set of equations of the plane-stress projected hyperelastoplastic model for large deformations is conveniently summarized in Box 2. The corresponding implicit stress integration algorithm is derived subsequently. 4. Implicit stress integration algorithm 4.1. Geometric pre-processor: trial elastic predictor The algorithm solves the stresses and strains at each Prandtl device, i.e. it is simply a repetition of the same algorithm for each device i = 1, ..., N . Assume that the total deformation gradient tensor for the plane stress problem t+∆t0 F is given in matrix notation at step t + ∆t (obtained from tentative displacements and passed to the stress-point material subroutine by the element subroutine). From previous computations, we also know the plastic deformation gradient t0 Fip , and the internal elastic strain and stress |e tensors t0 Eie , t Ti,α for each device i = 1, ..., N . We first consider the trial state, in which the 19

Box 2 Plane-stress constrained hyperelastoplastic model for large deformations 1: Elastic logarithmic strain tensor obtained from the multiplicative decomposition −1 T Fie = F Fip ⇒ Eie = 12 ln(Fie Fie ) 2: Elastic logarithmic strain rate tensor split E˙ ie = E˙ ie |F˙ p =0 + E˙ ie |E=0 = tr E˙ ie + ct E˙ ie ˙ i 3: In-plane generalized Kirchhoff stresses derived from the stored energy e P PN ∂Ei,α |e Tα = T∞,α + N i=1 Ti,α = T∞,α + i=1 Ti,α : ∂Eα 4: Evolution of the plastic flow at each Prandtl device i = 1, ..., N ; deviatoric and in the plane stress subspace ct e |e E˙ i,α = − 32 kγ˙ ii Ti,α 5: Yield function at each Prandtl device i = 1, ..., N already constrained in the plane-stress subspace |e |e |e fi (Ti,α , ki ) = Ti,α : Ti,α − 23 ki2 6: Loading/unloading conditions at the Prandtl devices |e |e γ˙ i ≥ 0, fi (Ti,α , ki ) ≤ 0, γ˙i fi (Ti,α , ki ) = 0 7: Relation totaland trial elastic logarithmic strains during plastic flow  between |ein-plane ¯ |e :T |e ) T ⊗( A i,α ˙ e = IS − |e i,α|e i,α E : tr E˙ e |e i,α

8:

¯ :T + 4 k2 k0 Ti,α :A i,α i,α 9 i i

i,α

In-plane internal continuum elastoplastic tangent matrix in the intermediate and reference configurations |e

¯ |ep = A ¯ |e − A i.α i,α ¯ i,α A

|e

|e

|e

¯ : T ) ⊗ (A ¯ :T ) (A i,α i,α i,α i,α

|e ¯ |e : T |e + 4 k 2 k 0 Ti,α : A i,α i,α 9 i i T e e ∂Ei,α ∂Ei,α ¯ |ep : = :A i,α ∂Eα F˙ p =0 ∂Eα F˙ p =0 i

i

20

plastic flow is frozen. In other words, we consider a purely elastic trial step defined by the following formulas: tr tr tr

Fip := t0 Fip

(83) p−1

Fie := t+∆t0 F t0 Fi

e Ei,α :=

tr |e Ti,α tr ki tr

:= :=

e → tr Fi,α

(84)

tr tr e 1 2 ln( Fi,α Fi,α ) t |e ¯ |e : (tr E e Ti,α + tr A i,α i,α t ki

(85)

eT

|e

e − t0 Ei,α )

|e

(86) (87)

|e

fi := f (tr Ti,α , tr ki ) = tr Ti,α : tr Ti,α −

2 tr 2 3 ki

(88)

Note that, since the third direction is a principal deformation direction, it is uncoupled, i.e. tr e  h i −1 Fi,α 0 tr e [ Fi ] := ≡ t+∆t0 F t0 Fip tr e 0 Fi,3 −1  t p t p t+∆t t+∆t 0 0 0 F11 0 F12 0 F11 0 F12 p p 0  = t+∆t0 F21 t+∆t0 F22 0  t0 F21 t0 F22 t p t+∆t 0 0 0 0 0 F33 0 F33    t+∆t t+∆t t+∆t t+∆t t p t p t p t p 1 0 F11 0 F12 0 F12 0 F21 0 F12 0 F11 − 0 F11 0 F22 − 0   t+∆t t+∆t t+∆t t+∆t p t p t p t p t p  det(t0 Fi,α )  0 F21 0 F12 0 F21 0 F22 − 0 F22 0 F21 0 F22 0 F11 − =  t+∆t  0 F33  0 t p 0 F33 (89) In the case tr fi ≤ 0, the step is elastic and the trial state is the solution at time t + ∆t. In the case tr fi > 0 , the trial elastic solution is no admissible for device i and the plastic correction must be computed for that device. 4.2. Plastic corrector: local Newton iterations If tr fi > 0, we consider a corrector step in which the external deformation is frozen and by internal evolution, the trial stress returns to the yield surface, i.e. |e

f ( t+∆t Ti,α ,

t+∆t

ki ) = 0

(90)

with ∆γi > 0. A plain backward-Euler (Closest Point Projection) algorithm is used here for the numerical implementation (we do not use the exponential mapping). The solution varie and the consistency parameter t+∆t γ . Herein, ables are the in-plane elastic strain t+∆t Ei,α i the residuals of both the elastic strains and the yield function are used to obtain the solutions. The first residue can be obtained through Relation (37) by employing the plain backward-Euler method Z t+∆t ˙ i |e t+∆t e t+∆t e tr e 3γ ρi,α = 0 Ei,α − Ei,α + 2 k Ti,α dτ i t ∆γ |e i e e = t+∆t0 Ei,α − tr Ei,α + 32 t+∆t t+∆t Ti,α → 0 (91) ki 21

where the internal elastic stress at time t + ∆t may be obtained at the intermediate configuration integrating Eq. (70) as t+∆t

|e

|e

Ti,α = t Ti,α +

t+∆t ¯ |e Ai,α

e e : (t+∆t0 Ei,α − t0 Ei,α )

(92)

The second residue is the yield function, i.e. using the form in Eq. (65) t+∆t

|e

|e

fi = t+∆t Ti,α : t+∆t Ti,α −

2 t+∆t 2 ki 3

→0

(93)

The two residues can be written in vector form as t+∆t e  t+∆t e  ρi,α 0 Ei,α Ri (E) = t+∆t → 0 with Ei = t+∆t γ fi i

(94)

The residual vector equation is solved using the Newton-Raphson method, where the solution is updated at iteration (j + 1) from the known values at iteration (j) by (j+1)

Ei

(j)

(j)

(j)

= Ei − [∇Ri ]−1 Ri

(95)

until (j+1)

||Ri

|| ≤ tol

(96)

where (•)(j) indicates quantities for iteration (j) at time step t + ∆t. For the first iteration, (0) e , t γ ]T . The Jacobian of the residual vector respect to the variables is we take Ei = [tr Ei,α i  ∂ρe (j)

∇Ri

=

i,α  ∂E e i,α   ∂fi e ∂Ei,α

∂ρei,α (j) ∂γi   ∂fi  ∂γi t+∆t (j) t+∆t k − ∆γ t+∆t k 0 |e i i i t+∆t 3 Ti,α  2 t+∆t k 2  i 0 − 43 t+∆t ki t+∆t ki

 3 ∆γi t+∆t ¯ |e S Ai,α ¯I + 2 t+∆t ki = |e ¯ |e 2 t+∆t Ti,α : t+∆t A i,α

(97)

e is updated at each iteration using the integration of The out-of-plane internal strain Ei,3 e Eq. (51) as follows—we consider here only the normal strain Ei,3

t+∆t e,(j+1) 0 Ei,3

=

t+∆t e,(j) 0 Ei,3

−

∂ 2 Ψei e ∂E e ∂Ei,3 i,3

!−1

(j) ∂ 2 Ψei : e ∂E e ∂Ei,3 i,α

e,(j+1)

: ( t+∆t0 Ei,α

−

t+∆t e,(j) 0 Ei,α )

t+∆t

(98) from which it is also immediate to determine the component of the elastic part of the defore = exp(E e ), whereas the update of the plastic one (from the corrector mation gradient Fi,3 i,3 elastic strains) follows from the incompressibility condition of plastic flow.

22

4.3. Geometric post-processor: update variables Once the the local iterations have converged,

t+∆t e t+∆t e t+∆t γ and i 0 Ei,α , 0 Ei,3 , t+∆t |e internal stress Ti,α in the updated

∆γi are

known. Then, we can obtain the in-plane intermediate configuration through Eq. (92), and we can perform a pull-back operation from the updated intermediate configuration to the reference configuration employing the proper mapping tensors. In order to obtain an algorithm with a simple linearization, we hinge on the the trial elastic configuration due to the fact that this trial configuration is frozen during the local plastic integration process. The stress in this trial configuration is defined as e , Ee ) e e , Ee ) dΨei (Ei,α ∂ t+∆t0 Ei,α d¯Ψei (Ei,α i,3 i,3 t+∆t |tr = : Ti,α = e e e ˙p ˙p d¯tr Ei,α dEi,α ∂ tr Ei,α Fi =0

'

t+∆t

|e Ti,α

¯S

:I =

t+∆t

t+∆t

Fi =0

|e Ti,α

(99)

where by d¯(•)/d¯(◦) we mean total derivative in the sense of scheme Eq. (48), but partial derivative in the sense that we still keep frozen the evolution of the configuration F˙ip = 0— i.e. a Lie-type derivative. An approximation has been made in this equation in the case of large nonproportional steps. No approximation takes place in the case of proportional steps (in such case logarithmic strains are additive) or small steps. Since we consider that in any case, for accuracy reasons, small plastic steps are needed in the non-proportional case, we can accept that the trial and final configurations for logarithmic strains are close enough, and e /∂ tr E e is approximately the 2D symmetric therefore, the in-plane mapping ∂ t+∆t0 Ei,α i,α p F˙i =0

identity tensor ¯IS . Taking again into account work-conjugacy and the fact that the third dimension is principal and uncoupled dE i,α Ti,α : E˙ i,α = Si,α : A˙ i,α ⇒ Si,α = Ti,α : dAi,α dΨ(Ai,α , Ai,3 (Ai,α )) dΨ(Ei,α , Ei,3 (Ei,α )) dE i,α ⇔ = : dAi,α dEi,α dAi,α

(100) (101)

Then, using the definition in Eq. (48) it is straightforward to perform another geometric mapping to obtain the respective second PK stresses in the trial configuration as follows t+∆t

|tr

Si,α =

e e ˙e d tr Ei,α d tr Ei,α dΨei |tr t+∆t |tr ¯ Ei,α : = T : = t+∆t Ti,α : M i,α ˙e tr e e e tr tr A d Ei,α d Ai,α i,α d Ai,α

(102)

˙e

e /d tr Ae ¯ Ei,α where M := d tr Ei,α i,α is a 2D fourth-order geometric mapping tensor obtained A˙ e i,α

directly from the spectral decomposition of the corresponding deformations (see Appendix A). The second PK internal stresses in the reference configuration are subsequently obtained mapping those from the trial configuration to the reference configuration as t+∆t

Si,α =

t+∆t

|tr Si,α

23

:

d tr Aei,α d t+∆t0 Aα

(103)

−T

−T

p p where d tr Aei,α /d t+∆t0 Aα = Fi,α Fi,α has been determined from the relationship in Eq. (35) using the in-plane box as explained in Eq. (89). The plastic deformation gradient tensor t+∆t Fip at step t + ∆t must be updated for the computation of the next incremental loading step. This computation does not present any difference with the 3D version because the third direction is always principal and rotationless. To obtain it, we need to update the elastic deformation gradient tensor t+∆t0 Fie . Here we assume that the elastic rotation tensor Re is not modified during the plastic flow, so we have t+∆t e t+∆t e tr e 0 Ri = Ri , and therefore 0 Fi can be updated using the right polar decomposition t+∆t e 0 Fi

= t+∆t0 Rie

t+∆t e 0 Ui

(104)

where the right stretch tensor t+∆t0 Uie can be obtained directly from the logarithmic elastic strain tensor as t+∆t0 Uie = exp( t+∆t0 Eie ). The plastic deformation gradient matrix is recovered as −1 t+∆t p Fi = t+∆t Fie t+∆t F (105) 4.4. Consistent tangent moduli For large deformations, if we follow the total Lagrangean scheme, we need to pass, at each global iteration for time step t + ∆t, the algorithmic elastoplastic tangent moduli t+∆t C = d t+∆ S/d t+∆t0 A to the element subroutines in order to preserve asymptotic quadratic convergence during the global equilibrium iterations. These moduli can be calculated by a pull-back operation from the trial intermediate configuration to the reference configuration. It is noted that the trial intermediate configuration of each device i remains fixed during the correction step for that device, and that the moduli obtained in this configuration are the |ep algorithmic counterpart of the continuum elastoplastic tangent Ai,α given in Eq. (79), which is defined as |e

|tr

t+∆t ¯ |tr Ai,α

:=

d t+∆t Ti,α e d tr Ei,α

'

|e

=

d t+∆t Ti,α e d t+∆t0 Ei,α

:

d t+∆t Ti,α e d tr Ei,α

e e d t+∆t0 Ei,α d t+∆t0 Ei,α t+∆t ¯ |e = A : i,α e e d tr Ei,α d tr Ei,α

(106)

e /d tr E e that operates As shown in Eq. (106), we need to determine the relation d t+∆t0 Ei,α i,α when plastic flow is taking place in Prandtl device i. We know that the first residue given in Eq. (91) t+∆t ρei,α = 0 must vanish between global iterations at step t + ∆t so it brings the following consistency equation d t+∆t ρei,α =0 (107) e d tr Ei,α

which results in e d t+∆t0 Ei,α + e d tr Ei,α

|e

t+∆t

Γi

d t+∆t Ti,α e d t+∆t0 Ei,α

e d t+∆t0 Ei,α d t+∆t Γ ¯S t+∆t |e : + T ⊗ i,α e e −I =0 d tr Ei,α d tr Ei,α

24

(108)

where we defined for convenience

:= 32 ∆γi /ki (γi ). Thus, we can factor out ! t+∆t Γ d |e i (109) : ¯IS − t+∆t Ti,α ⊗ tr e d Ei,α

t+∆t Γ (γ ) i i

e d t+∆t0 Ei,α t+∆t ¯ −1 = Ki e tr d Ei,α

¯ i is known once the local iterations have converged at each global iteration as—see where K Eq. (97) ∂ t+∆t ρei,α ∂

t+∆t e 0 Ei,α

|e

≡

t+∆t T 3 ∆γi d i,α S ¯ Ki = I + t+∆t = ¯IS + t+∆t e 2 ki d E 0 i,α

t+∆t ¯

t+∆t

¯ |e Γi A i,α

(110)

e is still unknown. In order to determine From Eq. (109), only the tensor d t+∆t Γi /d tr Ei,α it, we use that t+∆t fi = 0 between global iterations, so |e

t+∆t T e d t+∆t0 Ei,α 1 d t+∆t fi t+∆t |e d 2 i,α 0= : = T : − i,α e e t+∆t tr tr e 2 d Ei,α d Ei,α 3 d 0 Ei,α

t+∆t

ki ki0

d∆γi e d tr Ei,α

and using Eq. (109)   2 k0 k 4 d t+∆t Γi |e |e |e |e −1 i i ¯ |e : K−1 ) ¯ :K ¯ :T + = (T : A Ti,α : A i,α i,α i i,α i i,α e 3 3 − 2Γi ki0 t+∆t d tr Ei,α t+∆t

(111)

(112)

we can factor-out the remaining term d t+∆t Γi e d tr Ei,α

2 0 ¯ i : T |e + 4 ki ki 0 :D i,α 3 3−2Γi k |e ¯i Ti,α : D

=

|e

Ti,α

i

(113)

t+∆t

¯i = A ¯ |e : K ¯ −1 is defined as the in-plane elastic algorithmic moduli tensor. It is noted where D i,α i ¯ i → ¯IS and D ¯i → A ¯ |e . Substituting this formula into that when ∆γi → 0 the tensor t+∆t K i,α

Eq. (109), we have e d t+∆t0 Ei,α e d tr Ei,α

 ¯ −1 : ¯IS − =K i

  2 k0 k |e 4 i i ¯ : Di : Ti,α + 3 3−2Γi k0

|e Ti,α |e

Ti,α

¯ i : T |e ) ⊗ (D i,α

i

(114) t+∆t

which clearly approaches the relation presented in Eq. (77) in the continuum limit, with Γi → 0. Finally, substituting this last equation into Eq. (106) yields   ¯ i : T |e ) ⊗ (D ¯ i : T |e ) ( D t+∆t ¯ |tr i,α i,α ¯i −  (115) Ai,α = D 2 0 |e ¯ i : T|e + 4 ki ki 0 Ti,α : D i,α 3 3−2Γi k i

t+∆t

¯ |tr → ¯ i → ¯IS and D ¯i → A ¯ |e , so t+∆t A where it is immediate to check that when Γi → 0 then K i,α i,α t+∆t ¯ |ep Ai,α , which is the continuum tangent given in Eq. (79). The internal elastoplastic 25

consistent tangent in the trial configuration for quadratic strains can be defined as ! |tr e d t+∆t Si,α d tr Ei,α d t+∆t ¯ |tr t+∆t |tr Ci,α = = tr e Ti,α : tr e d tr Aei,α d Ai,α d Ai,α !T e e e d tr Ei,α d tr Ei,α d2 tr Ei,α t+∆t ¯ |tr t+∆t |tr = : A : + T : i,α i,α d tr Aei,α d tr Aei,α d tr Aei,α d tr Aei,α

(116)

where the sixth-order tensor is given, for example, in Ref. [68]. The final pull-back operation t+∆t ¯ |tr t+∆t ¯ over Ci to give Ci in the reference configuration becomes straightforward by a mapping between the fixed trial configuration and reference configuration as !T |tr d t+∆t Si,α d tr Aei,α d tr Aei,α d t+∆t Si,α = Ci = t+∆t : : t+∆t d tr Aei,α d d t+∆t0 Ai,α d 0 Ai,α 0 Ai,α !T d tr Aei,α d tr Aei,α t+∆t ¯ |tr C : : = i,α d t+∆t0 Ai,α d t+∆t0 Ai,α

t+∆t ¯

(117)

The kinematic tangent moduli matrix in the reference configuration can be obtained as !T

d t+∆t0 Eα t+∆t d2 t+∆t0 Eαe + T : ∞ d t+∆t0 Aα d t+∆t0 Aα d t+∆t0 Aα (118) and, finally, the consistent in-plane elastoplastic tangent matrix in the reference configuration, for global equilibrium iterations during finite element analyses, is t+∆t ¯

C∞

d t+∆t S∞,α = t+∆t = d 0 Aα

d t+∆t0 Eα d t+∆t0 Aα

t+∆t ¯

:

C=

t+∆t ¯

A∞,α :

t+∆t ¯

C∞ +

N X t+∆t

¯i C

(119)

i=1

A summary of the implicit stress integration algorithm is given in Boxes 3 and 4. Remark : Although we present the formulation and the examples below for the case of isotropic elasticity and isotropic yield functions, the present approach may be applied to anisotropic hyperelasticity and anisotropic yield functions with very minor changes. For including anisotropic elasticity, no change is needed, singe we have already included the possibility of general hyperelasticity. For including anisotropic yield functions, we just need to change the flow rule to a form ct

|e e E˙ i,α = −γ˙ i NT : Ti,α

(120)

where NT is an anisotropy tensor (e.g. a Hill-type tensor). Then, changes in the algorithm are minimal, as in the 3D case; see details in Ref. [60].

26

Box 3 Plane stress-projected implicit integration algorithm 1: Trial elastic predictor: |e Given t+∆t0 F , t0 Fip , t0 Eie , t Ti,α and t0 γi , compute the trial elastic predictor state for each device i T tr p Fi := t0 Fip ; tr Fie := t+∆t0 F t0 Fip ; tr Eie := 21 ln( tr Fie tr Fie ); tr T |e i,α

|e

|e

e − t Ee ) := t Ti,α + tr Ai,α : ( tr Ei,α 0 i,α 2: Check yield conditions and test for plastic loading at each device i tr ki := t ki |e |e |e tr fi := fi ( tr Ti,α , tr ki ) = tr Ti,α : tr Ti,α − 32 tr ki2 IF tr fi ≤ 0 THEN Elastic step: set t+∆t (•) = tr (•) and proceed to next device ELSE Plastic step: Proceed to step 3 ENDIF 3: Return mapping for device i

(i) Local Newton  iterations:plastic corrector t+∆t ρ i,α Solve Ri (E) = =0 t+∆t f i e,(0)

(0)

(0)

e , ∆γ with initial values t+∆t Ei,α = tr Ei,α = 0, and t+∆t ki i (ii) Update the variables t+∆t T |tr ' t+∆t T |e = t T |e + t+∆t A ¯ |e : ( t+∆t E e − t E e ) 0 i,α 0 i,α i,α i,α i,α i,α tr E e d i,α |tr |tr t+∆t S t+∆t T i,α = i,α : tr e d Ai,α dtr Aei,α |tr t+∆t S t+∆t Si,α : t+∆t i,α = d Ai,α t+∆t F p = t+∆t F e−1 t+∆t F i i

= ki ( t γi )

P (iii) After all devices, add all contributions: t+∆t Sα = t+∆t S∞,α + N i=1 (iv) Compute the in-plane elastoplastic tangent matrix — see Box 4

27

t+∆t S |e i,α

Box 4 In-plane consistent elastoplastic tangent 1: For each device, compute the internal elastoplastic tangent in the trial configuration for logarithmic  strains  |e |e ¯ ¯ ( D :T )⊗( D :T ) i i,α i i,α t+∆t A ¯ |tr = D ¯i −  2 k0 i,α k 4 |e ¯ |e i i Ti,α :D i :Ti,α + 3 3−2Γ k 0 i i

t+∆t A ¯ |tr i,α

Map ration. 3: Pull back 2:

t+∆t ¯

Ci =

4:

to quadratic measures

t+∆t t+∆t C ¯ |tr i

|tr

= d t+∆t Si,α /d tr Aei,α in the trial configu-

t+∆t ¯ |tr Ci to the reference configuration d tr Aei,α t+∆t |tr d tr Aei,α ¯ : : C i d t+∆t0 Aα d t+∆t0 Aα

Compute the consistent elastoplastic tangent matrix adding all contributions t+∆t ¯

C=

with

t+∆t ¯

C∞ +

t+∆t ¯

C∞ =

PN

i=1

t+∆t ¯

Ci

d t+∆t Eα t+∆t ¯ d t+∆t Eα t+∆t d2 t+∆t0 Eαe : A : + T : ∞,α ∞ d t+∆t0 Aα d t+∆t0 Aα d t+∆t0 Aα d t+∆t0 Aα

5. Numerical examples In this section, we show some demonstrative examples to verify that the proposed planestress projected formulation and its numerical implementation are efficient and give similar results to the 3D model when the thickness is small respect to the in-plane dimensions, resulting therefore in congruent formulations. Furthermore, the classical approach for plane stress constraint at the Gauss point level based on a nested loop [2, 3] is also implemented (see Appendix B for a summary of the iterative algorithm, after which kinematic mappings are applied) in order to compare with the proposed plane-stress projected formulation. 5.1. Bending of an end-loaded cantilever The first problem is the plastic collapse of an end-loaded cantilever beam with dimensions shown in Fig. 2. The material for this problem is perfectly plastic with the following properties: Young modulus E = 210 GPa, Poisson coefficient ν = 0.3 and yield stress σY = 240 MPa. In this problem, only one Prandtl device is needed. The cantilever is modeled for both 2D and 3D analysis in order to verify the functionality of the proposed formulation. For the 2D analysis, the cantilever is discretized by 2 × 50 u/p fully integrated (3 × 3 Gauss integration) Q2/P1-9/3 finite elements [1] using the plane-stress projected formulation, while for the 3D analysis, the same discretization is used with the u/p fully integrated (3 × 3 × 3 Gauss integration) Q2/P1–27/4 brick finite elements without the plane stress constraint (3D model). Note that for the plane stress model we could use a standard quadratic finite element formulation since no volumetric locking effect is expected and quadratic elements do not lock in bending. The numerical solution for both 2D and 3D analyses is obtained by applying a vertical prescribed displacement of 6 mm at the loaded node, and the equivalent applied force can be determined as the total reaction force at any side of the beam. It is noted 28

that the vertical deflection is 12 % respect to the beam length, which means a significant large displacement, that further results in large plastic strains near the support of the beam. The accumulated plastic deformation is shown in Fig. 3 for all analyses. There is a small difference in the plastic zone in the sense that in the 3D analysis the plastic deformations are slightly more localized near the support due to the tridimesionality effect. Note that both element types use also a different number of pressure interpolation variables. The plastic zone obtained in 2D analysis by the proposed plane-stress projected model and the nested loop model are similar. There is more plastic deformation with the nested loop model. In general, despite the expected small differences, congruent results can be observed.

Figure 2: End-loaded cantilever: geometric sketch and FE model (all dimensions in mm)

The analytical limit load for triggering plasticity is given by Flim =

σY bh2 4L

(121)

where b, h are the width and height of the rectangular cross beam section, respectively, and L is the length of the cantilever. The above solution is derived for the infinitesimal framework in Section 4.4.2 of Ref. [74]. For the present geometry and material properties, the analytical limit load is F = 30 kN. The equivalent force applied at the end of the cantilever is obtained for all analyses and plotted versus the vertical deflection at the loaded node are shown in Fig. 4. It can be observed that 2D analysis with the proposed plane-stress projected model and 3D analysis give similar results and that the plastic load is also in agreement with the analytical one, even though in our calculation the finite strain framework is used; note that large displacements and also large plastic strains are reached as shown in Fig. 3. The result obtained by the nested loop model gives slightly different values than the other two analyses as expected, since from Fig. 3 the nested model results in more plastic deformation than the other two. It can be observed that the values given by the nested loop algorithm are on the softer side during hardening respect to the 3D solution, whereas the solution by the 2D plane stress projected algorithm is on the stiffer side respect to the 3D solution. Noteworthy, this latter solution would be the one expected in the 2D analysis, because in this latter case the plane stress constrain is fully enforced, whereas in the 3D case the out-of-plane stresses are not fully constrained. This difference may be due to the fact that all the addends of the stresses derived from our plane stress constrained algorithm are in-plane by construction, (as 29

Figure 3: End-loaded cantilever. Accumulated plastic deformation at the vertical prescribed deflection of 8 mm: (a) results obtained in 2D analysis with proposed plane-stress projected model, (b)results obtained in 2D analysis with nested loop model, (c) results obtained in 3D analysis.

obtained from the dissipation equation), whereas the nested loop algorithm only constrains the total, final stress, not the individual addends. This is consistently observed in all the examples. Typical convergence rates during the global Newton-Raphson iterations for step 100 are shown in Fig. 5, where it is seen that the convergence rate was a bit faster in the proposed plane-stress projected model. Furthermore, it is noted that the the total time consumed in the 2D analysis with the plane-stress projected formulation is 1434s, lower than the 2D nested loop model of 1913s and is much efficiently than the 3D model of 6723s. 5.2. Extension of a strip with a circular hole We perform in this example a numerical simulation of the extension of a strip of thickness 1 mm with a central circular hole, a benchmark problem addressed in many works, e.g. [14, 3] and therein references. The structure is subjected to a loading simulated via imposed displacement in the vertical direction, see Fig. 6. Only a quarter part of the plate is modeled using the u/p fully integrated (3 × 3 Gauss integration) Q2/P1-9/3 finite elements for the 2D analysis because of symmetry considerations; while the mixed u/p fully integrated (3 × 3 × 3 Gauss integration) Q2/P1 - 27/4 brick finite elements are used for the 3D analysis in order to compare with the results obtained with the plane-stress projected algorithm. The material used for this simulation is a DP600 dual-phase steel with nonlinear kinematic hardening and its parameters obtained from the experimental curve are presented in Ref. [36]. We herein discretized the uniaxial behavior in 11 internal Prandtl devices which properties are 30

40

Applied load (kN)

30

20

Force Force Force Anal.

10

0

0

1

in 2D plane-stress projected model in 2D nested loop model in 3D limit load in small strain

2

3

4

5

6

Vertical deflection at loaded node (mm) Figure 4: End-loaded cantilever: Comparison of the load-deflection diagram using the 3D formulation [63] and the present 2D constrained formulation. Dash line corresponds to the limit load for the infinitesimal strains theory

2D plane-stress projected model

2D nested loop model 104

104

10−4

101

enorm

rnorm

3D model

10−12

10−2

10−5

10−20 1

2

3

4

5

1

Number of global iterations

2

3

4

5

Number of global iterations

Figure 5: Bending of an end-loaded cantilever: Convergence rate for global (equilibrium) Newton-Raphson iterations (rnom: residual in force; enorm: residual in energy) for step 100

computed following the scheme described in Ref. [63]. The material parameters are listed in Table 1. The prescribed displacement is up to a total of 0.4 mm, which corresponds to an average deformation of 2.5% (substantially larger deformations are obviously obtained in the boundaries of the hole). The solution has been obtained with our in-house finite element code Dulcinea. The evolution of the accumulated plastic (logarithmic) strain at different stages of the loading process is shown in Fig. 7 for both 2D analyses (with the proposed plane-stress projected

31

Cauchy stress σ11 (MPa)

800

600

400

200 Exp. data 0

0

5 · 10−2

0.1

0.15

0.2

Logarithmic strain E11 (mm) (a)

(b)

Figure 6: Strip with circular hole: (a) Geometry and FE model (all dimensions in mm); (b) Cauchy axial stress versus logarithmic axial strain for DP600 steel in Ref. [36]

formulation and with the nested loop - based approach) and for the 3D analysis. It can be observed that the evolution of the accumulated plastic strain obtained with the plane-stress projected model is very close to the one obtained with the nested loop - based approach and the 3D model at each stage. Indeed, the maximum accumulated plastic strain reached at the final prescribed displacement is similar: 0.2131 for 2D plane stress-projected model, 0.2308 for 2D nested loop - based model and 0.2225 for the 3D model. A difference that can be expected due to the 3D effect for the given thickness in the example (plane stress conditions are not physically exact) and the different pressure interpolation degrees of freedom needed for passing the inf-sup condition [1]. The von Mises stress bandplots are shown in Fig. 8. Similar conclusions are obtained from this figure, showing the agreement between both algorithms. It is interesting to note also in this example that the solution for the 2D nested loop model is in the flexible side, whereas the one using the 2D plane stress projected formulation is in the stiff side. As mentioned, since the 2D problem fully constrains the solution, the latter result is to be expected. To investigate which solution would be the limit as the thickness goes to zero, we have repeated the 3D analysis with 1/2 and 1/4 of the original thickness. In Table 2 we show the maximum accumulated plastic strains. It is seen that as the thickness goes to zero, the solution gets closer to the proposed 2D algorithm and farther from the nested loop algorithm. In Fig. 9 we show the forces and energy convergence rates for step 100. These rates show that, as the 3D formulation, the plane-stress algorithm also presents the asymptotic quadratic convergence rates of Newton procedures. Finally, we note that the 2D analysis (total time consumed of 6344s) employing the herein proposed plane-stress projected formulation is four times faster (and needs substantially less memory) than the 3D analysis (total time consumed of 23061s) and is 1.5 times faster than the 2D analysis with the nested loop- based approach

32

Table 1: Material parameters of DP600 steel for proposed material model

Young modulus E = 206 GPa Poisson coef. ν = 0.28 External shear moduli 2µ∞ = 439.02 MPa Internal Prandtl devices parameters 2µ1 = 102746.01 MPa No. 1 k1 = 308.24 MPa 2µ2 = 1200.65 MPa No. 2 k2 = 17.35 MPa 2µ3 = 513.53 MPa No. 3 k3 = 15.67 MPa 2µ4 = 663.15 MPa No. 4 k4 = 28.98 MPa 2µ5 = 720.06 MPa No. 5 k5 = 41.82 MPa 2µ6 = 518.73 MPa No. 6 k6 = 39.97 MPa 2µ7 = 382.00 MPa No. 7 k7 = 42.40 MPa 2µ8 = 549.39 MPa No. 8 k8 = 70.79 MPa 2µ9 = 296.70 MPa No. 9 k9 = 50.20 MPa 2µ10 = 32.95 MPa No. 10 k10 = 6.11 MPa 2µ11 = 21.97 MPa No. 11 k11 = 4.83 MPa Table 2: Convergence of the 3D model towards the 2D plane stress model

Model 3D with thickness of 1 mm 3D with thickness of 0.5 mm 3D with thickness of 0.25 mm 2D plane stress projected algorithm (proposal) 2D plane stress with nested loop algorithm

Max. plastic strain 0.2225 0.2208 0.2185 0.2131 0.2308

model (total time consumed of 9739s). Note that these computational times include all tasks in the program, most of which are identical to both 2D analyses, so the relative differences between both 2D formulations are larger.

33

Figure 7: Evolution of the accumulated plastic (logarithmic) strains at different stages of the loading process: (a) results obtained with 2D plane-stress projected algorithm; (b) results obtained with the 2D nested loop algorithm; (c) results obtained with full 3D model.

34

Figure 8: The von Mises stress at different stages of the loading process: (a) Results obtained with 2D plane-stress model; (b) Results obtained with 2D nested loop approach; (c) Results obtained with full 3D model.

35

2D plane-stress projected model

2D nested loop model 101

100

10−6

rnorm

enorm

103

3D model

10−3

10−13

10−6 1

2

3

4

10−20

5

Number of global iterations

1

2

3

4

5

Number of global iterations

Figure 9: Extension of strip with circular hole: Convergence rate for global (equilibrium) Newton-Raphson iterations (rnom: residual in force; enorm: residual in energy) for step 100

5.3. Biaxial cyclic loading of a strip with a circular hole In this example, we consider the same geometry and material as the previous one but subjected to a biaxial cyclic loading history shown in Fig. 10. As in the previous example, nonlinear kinematic hardening is considered without isotropic hardening as to test the preservation of Masing’s rules under multiaxial, nonproportional loading even when using the proposed plane-stress projected model. The numerical calculations are performed again using both the 3D and the 2D plane-stress formulations. The evolution of the von Mises stress at different states of the cyclic loading process is shown in Fig. 11. In Figure 10, several loading instants are found to be equivalent regarding Masing effect. For times t = 0.2s and t = 1s, the loading at both sides is the same and the trajectory from previous steps is also the same; hence very similar stress bandplots are to be expected for a formulation that preserves Masing’s rules. It is seen that indeed, there is no relevant difference between both instants in Fig. 11. A similar comparison may be done between instants t = 0.4s and t = 1.2s. Hence, the plane-stress projected formulation with nonlinear kinematic hardening is capable of fulfilling Masing’s behavior in a multiaxial, nonproportional setting. Of course, additional isotropic hardening may be included as described in the previous sections. Finally, the accumulated plastic strain is shown in Fig. 12. 6. Conclusions In this paper, a plane-stress constrained formulation for large strain multiplicative hyperelastoplasticity which includes nonlinear kinematic hardening and is derived from the principle of the maximum dissipation, is presented. The formulation is based on the elastic logarithmic strain rate corrector concept and results in continuum and algorithmic frameworks which are parallel to both the small strains setting and to the full 3D version. The consistent tangent also preserves the same simplicity and parallelism, and gives second order convergence rates 36

Prescribed displacement (mm)

Upper side

Left side

0.4

0.8

0.4

0.2

0

−0.2

−0.4

0

0.2

0.6

1

1.2

Time (s) Figure 10: Biaxial cyclic loading functions for strip with circular hole

in Newton global equilibrium iterations. Furthermore, if nonlinear purely kinematic hardening is considered, Masing’s behavior is also replicated in the plane stress-projected model, as an inherited property from the 3D version. Noteworthy, the presented work is the first plane-stress projected formulation (although we prefer in this case to refer to it as a plane-stress constrained formulation) for multiplicative hyperelastoplasticity which considers nonlinear kinematic hardening. The soundness, robustness and remarkable simplicity of the formulation and its algorithmic implementation, within the context of multiplicative plasticity and nonlinear kinematic hardening, further demonstrates the usefullness of the recently proposed framework in terms of elastic logarithmic corrector rates teamed with Besseling rehological models. Acknowledgement Partial financial support for this work has been given by Agencia Estatal de Investigaci´ on of Spain under grants DPI2015-69801-R and PGC2018-097257-B-C32. [1] K.-J. Bathe, Finite Element Procedures, 2nd Ed, Klaus-J¨ urgen Bathe, 2014. [2] A. Millard, Numerical Algorithms for Plane Stress Elastoplasticity: Review and Recommendation, in: D. R. J. Owen, E. O˜ nate (Eds.), 4th International conference, Computational plasticity: fundamentals and applications, Pineridge, Barcelona; Spain, 1995, pp. 237–248. [3] E. A. de Souza Neto, D. Peric, D. R. J. Owen, Computational Methods for Plasticity, John Wiley & Sons, Ltd, Chichester, UK, 2008. [4] D. M. Tracey, C. E. Freese, Adaptive load incrementation in elastic-plastic finite element analysis, Computers and Structures 13 (1-3) (1981) 45–53.

37

Figure 11: Biaxial cyclic loading. Evolution of von Mises stress at different stages of the cyclic loading process: (a) 2D analysis with plane-stress projected formulation , (b) 3D analysis

[5] R. H. Dodds, Numerical techniques for plasticity computations in finite element analysis, Computers and Structures 26 (5) (1987) 767–779. [6] R. G. Whirley, J. O. Hallquist, G. L. Goudreau, An assessment of numerical algorithms 38

Figure 12: Biaxial cyclic loading. Evolution of plastic zones at different stages of the cyclic loading process (proposed formulation).

for plane stress and shell elastoplasticity on supercomputers, Engineering Computations 6 (2) (1989) 116–126. [7] S. Klinkel, S. Govindjee, Using finite strain 3D-material models in beam and shell elements, Engineering Computations 19 (3-4) (2002) 254–271. [8] R. de Borst, The zero-normal-stress condition in plane-stress and shell elastoplasticity, Communications in Applied Numerical Methods 7 (1) (1991) 29–33. [9] P. Fuschi, D. Peri´c, D. R. Owen, Studies on generalized midpoint integration in rateindependent plasticity with reference to plane stress J2-flow theory, Computers and Structures 43 (6) (1992) 1117–1133. [10] P. B. Louren¸co, R. de Borst, J. G. Rots, A plane stress softening plasticity model for orthotropic materials, International Journal for Numerical Methods in Engineering 40 (21) (1997) 4033–4057. 39

[11] F. J. Mont´ ans, Implicit plane stress algorithm for multilayer J2-plasticity using the Prager-Ziegler translation rule, International Journal for Numerical Methods in Engineering 59 (3) (2004) 409–418. [12] N. Valoroso, L. Rosati, Consistent derivation of the constitutive algorithm for plane stress isotropic plasticity. Part I: Theoretical formulation, International Journal of Solids and Structures 46 (1) (2009) 74–91. [13] N. Valoroso, L. Rosati, Consistent derivation of the constitutive algorithm for plane stress isotropic plasticity. Part II: Computational issues, International Journal of Solids and Structures 46 (1) (2009) 92–124. [14] J. C. Simo, T. J. Hughes, Computational inelasticity, Vol. 7, Springer Science & Business Media, 2006. [15] J. C. Simo, R. L. Taylor, A return mapping algorithm for plane stress elastoplasticity, International Journal for Numerical Methods in Engineering 22 (3) (1986) 649–670. [16] P. Jetteur, Implicit integration algorithm for elastoplasticity in plane stress analysis, Engineering Computations 3 (1986) 251–253. [17] E. Ramm, A. Matzenmiller, Copmputational aspects of elasto-plasticity in shell analysis, in: O. E. Owen DRJ, Hinton E (Ed.), Computational Plasticity: Models, Software and Applications, 1987. [18] P. J. Armstrong, C. O. Frederick, A mathematical representation of the multiaxial Bauschinger effect, Tech. rep., Report RD/B Berkeley, (also available in Materials at High Temperatures 2007; 24: 1-26) (1966). [19] J.-L. Chaboche, Time-independent constitutive theories for cyclic plasticity, International Journal of Plasticity 2 (2) (1986) 149–188. [20] N. Ohno, M. Tsuda, T. Kamei, Elastoplastic implicit integration algorithm applicable to both plane stress and three-dimensional stress states, Finite Elements in Analysis and Design 66 (2013) 1–11. [21] O. Hopperstad, S. Remseth, A return mapping algorithm for a class of cyclic plasticity models, International Journal for Numerical Methods in Engineering 38 (1995) 549–564. [22] E. Kirchner, S. Reese, P. Wriggers, A finite element method for plane stress problems with large elastic and plastic deformations, Communications in Numerical Methods in Engineering 13 (1997) 963–976. [23] M. Kojic, K.-J. Bathe, Inelastic Analysis of Solids and Structures, Springer-Verlag Berlin Heidelberg, 2005. [24] P. Wriggers, Nonlinear finite element methods, Springer, 2008.

40

[25] R. de Borst, M. A. Crisfield, J. J. C. Remmers, C. V. Verhoosel, Non-Linear Finite Element Analysis of Solids and Structures, John Wiley & Sons, Ltd, Chichester, UK, 2012. [26] T. J. R. Hughes, J. Winget, Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis, International Journal for Numerical Methods in Engineering 15 (12) (1980) 1862–1867. [27] W. Rolph III, K. Bathe, On a large strain finite element formulation for elasto-plastic analysis, in: K. J. Willam (Ed.), Constitutive Equations, Macro and Computational Aspects. Winter Annual Meeting of ASME, ASME: New York, 1984, pp. 131–147. [28] G. C. Johnson, D. J. Bammann, A discussion of stress rates in finite deformation problems, International Journal of Solids and Structures 20 (8) (1984) 725–737. [29] A. Meyers, H. Xiao, O. Bruhns, Choice of objective rate in single parameter hypoelastic deformation cycles, Computers & Structures 84 (17-18) (2006) 1134–1140. [30] A. Shutov, J. Ihlemann, Analysis of some basic approaches to finite strain elastoplasticity in view of reference change, International Journal of Plasticity 63 (2014) 183– 197. [31] H. Xiao, O. Bruhns, A. Meyers, Hypo-elasticity model based upon the logarithmic stress rate, Journal of Elasticity 47 (1) (1997) 51–68. [32] B. Bernstein, Hypo-elasticity and elasticity, Archive for Rational Mechanics and Analysis 6 (1) (1960) 89–104. [33] B. Bernstein, Relations between hypo-elasticity and elasticity, Transactions of the Society of Rheology 4 (1) (1960) 23–28. [34] J. Simo, On the computational significance of the intermediate configuration and hyperelastic stress relations in finite deformation elastoplasticity, Mechanics of Materials 4 (3-4) (1985) 439–451. [35] J. Simo, M. Ortiz, A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations, Computer Methods in Applied Mechanics and Engineering 49 (2) (1985) 221–245. [36] I. N. Vladimirov, M. P. Pietryga, S. Reese, Prediction of springback in sheet forming by a new finite strain model with nonlinear kinematic and isotropic hardening, Journal of Materials Processing Technology 209 (8) (2009) 4062–4075. [37] I. N. Vladimirov, M. P. Pietryga, S. Reese, Anisotropic finite elastoplasticity with nonlinear kinematic and isotropic hardening and application to sheet metal forming, International Journal of Plasticity 26 (5) (2010) 659–687.

41

´ Caminero, F. J. Mont´ [38] M. A. ans, K. J. Bathe, Modeling large strain anisotropic elastoplasticity with logarithmic strain and stress measures, Computers and Structures 89 (1112) (2011) 826–843. [39] J. L¨oblein, J. Schr¨ oder, F. Gruttmann, Application of generalized measures to an orthotropic finite elasto-plasticity model, Computational Materials Science 28 (3-4) (2003) 696–703. [40] C. Miehe, A formulation of finite elastoplasticity based on dual co-and contra-variant eigenvector triads normalized with respect to a plastic metric, Computer Methods in Applied Mechanics and Engineering 159 (3-4) (1998) 223–260. [41] C. Miehe, N. Apel, M. Lambrecht, Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials, Computer Methods in Applied Mechanics and Engineering 191 (47-48) (2002) 5383–5425. [42] P. Papadopoulos, J. Lu, A general framework for the numerical solution of problems in finite elasto-plasticity, Computer Methods in Applied Mechanics and Engineering 159 (1-2) (1998) 1–18. [43] P. Papadopoulos, J. Lu, On the formulation and numerical solution of problems in anisotropic finite plasticity, Computer Methods in Applied Mechanics and Engineering 190 (37-38) (2001) 4889–4910. [44] C. Sansour, W. Wagner, Viscoplasticity based on additive decomposition of logarithmic strain and unified constitutive equations: Theoretical and computational considerations with reference to shell applications, Computers & Structures 81 (15) (2003) 1583–1594. [45] F. Aldakheel, P. Wriggers, C. Miehe, A modified Gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling, Computational Mechanics 62 (4) (2018) 815–833. [46] P. Neff, I.-D. Ghiba, Loss of ellipticity for non-coaxial plastic deformations in additive logarithmic finite strain plasticity, International Journal of Non-Linear Mechanics 81 (2016) 122–128. [47] T. Brepols, I. N. Vladimirov, S. Reese, Application and evaluation of hyper-and hypoelastic-based plasticity models in the FE simulation of metal forming processes, Proceedings in Applied Mathematics and Mechanics (PAMM) 13 (1) (2013) 153–154. [48] T. Brepols, I. N. Vladimirov, S. Reese, Numerical comparison of isotropic hypo- and hyperelastic-based plasticity models with application to industrial forming processes, International Journal of Plasticity 63 (2014) 18–48. [49] J. C. Simo, A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition: Part i. continuum formulation, Computer Methods in Applied Mechanics and Engineering 66 (2) (1988) 199–219. 42

[50] J. C. Simo, A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part ii: computational aspects, Computer Methods in Applied Mechanics and Engineering 68 (1) (1988) 1–31. [51] J. C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, Computer Methods in Applied Mechanics and Engineering 99 (1) (1992) 61–112. [52] G. Weber, L. Anand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids, Computer Methods in Applied Mechanics and Engineering 79 (2) (1990) 173–202. [53] A. L. Eterovic, K.-J. Bathe, A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures, International Journal for Numerical Methods in Engineering 30 (6) (1990) 1099–1114. [54] A. Cuiti˜ no, M. Ortiz, A material-independent method for extending stress update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics, Engineering Computations 9 (4) (1992) 437–451. [55] A. Lion, Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models, International Journal of Plasticity 16 (5) (2000) 469–494. [56] D. L. Henann, L. Anand, A large deformation theory for rate-dependent elastic–plastic materials with combined isotropic and kinematic hardening, International Journal of Plasticity 25 (10) (2009) 1833–1878. [57] M. E. Gurtin, E. Fried, L. Anand, The mechanics and thermodynamics of continua, Cambridge University Press, 2010. [58] A. Shutov, R. Kreißig, Finite strain viscoplasticity with nonlinear kinematic hardening: Phenomenological modeling and time integration, Computer Methods in Applied Mechanics and Engineering 197 (21-24) (2008) 2015–2029. [59] M. Latorre, F. J. Mont´ ans, A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate, Applied Mathematical Modelling 55 (2018) 716–740. [60] M. Sanz, F. J. Mont´ ans, M. Latorre, Computational anisotropic hardening multiplicative elastoplasticity based on the corrector elastic logarithmic strain rate, Computer Methods in Applied Mechanics and Engineering 320 (2017) 82–121. [61] E. Kr¨ oner, Allgemeine kontinuumstheorie der versetzungen und eigenspannungen, Archive for Rational Mechanics and Analysis 4 (1) (1959) 273. [62] E. H. Lee, Elastic-plastic deformation at finite strains, Journal of Applied Mechanics 36 (1) (1969) 1–6. 43

[63] M. Zhang, F. J. Mont´ ans, A simple formulation for large-strain cyclic hyperelastoplasticity using elastic correctors. theory and algorithmic implementation, International Journal of Plasticity 113 (2019) 185–217. [64] M. L. Wilkins, Calculation of elastic-plastic flow, Tech. rep., Report, California Univ Livermore Radiation Lab (1963). [65] J. C. Simo, Numerical analysis and simulation of plasticity, Vol. 6, Elsevier, 1998, pp. 183–499. [66] M. Mi˜ nano, M. A. Caminero, F. J. Mont´ans, On the numerical implementation of the closest point projection algorithm in anisotropic elasto-plasticity with nonlinear mixed hardening, Finite Elements in Analysis and Design 121 (2016) 1–17. [67] R. I. Borja, Plasticity. Modeling & Computation., Springer Heidelberg, New York, 2013. [68] M. Latorre, F. J. Mont´ ans, Stress and strain mapping tensors and general workconjugacy in large strain continuum mechanics, Applied Mathematical Modelling 40 (56) (2016) 3938–3950. [69] J. Crespo, M. Latorre, F. J. Mont´ans, WYPIWYG hyperelasticity for isotropic, compressible materials, Computational Mechanics 59 (1) (2017) 73–92. [70] J. Crespo, F. J. Mont´ ans, A continuum approach for the large strain finite element analysis of auxetic materials, International Journal of Mechanical Sciences 135 (2018) 441–457. [71] M. Latorre, F. J. Mont´ ans, Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials, European Journal of Mechanics-A/Solids 53 (2015) 99–106. [72] E. De Rosa, M. Latorre, F. J. Mont´ans, Capturing anisotropic constitutive models with wypiwyg hyperelasticity: and on consistency with the infinitesimal theory at all deformation levels, International Journal of Non-Linear Mechanics 96 (2017) 75–92. [73] M. Latorre, F. J. Mont´ ans, Fully anisotropic finite strain viscoelasticity based on a reverse multiplicative decomposition and logarithmic strains, Computers & Structures 163 (2016) 56–70. [74] J. Lubliner, Plasticity theory, Dover Publications, New York, 2008.

44

Appendix A. In-plane mapping tensors For the plane stress problem, the spin tensor of the material principal directions can be defined as: Ω=

3 X X i=1 j6=i

=

2 X X

Ωij ni ⊗ nj =

2 X X

Ωij ni ⊗ nj +

i=1 j6=i

2 X i=1

2 X 0 *n0 ⊗ n +  >n ⊗ n  Ω Ω i3 i 3 3j j  3 j=1

Ωij ni ⊗ nj = Ω12 (n1 ⊗ n2 − n2 ⊗ n1 )

(A.1)

i=1 j6=i

that depends only on the in-plane terms, and has an in-plane representation in principal ¯ i as in-plane directions n   0 −Ω 12 ¯ = Ω12 [n ¯1 ⊗ n ¯2 − n ¯2 ⊗ n ¯ 1] = [Ω] (A.2) Ω12 0 Therefore, following the mathematical development in Ref. [68], the tensor that maps the in-plane logarithmic strain Eα to the in-plane Green-Lagrange strain Aα can be determined as: ¯ E˙ α := dEα = 1 M ¯ 1 ⊗ M1 + 1 M ¯2 ⊗ M ¯ 2 + 4(log λ2 − log λ1 ) M ¯S ⊗M ¯S M 12 12 2 2 A˙ α dAα λ1 λ2 λ22 − λ21

(A.3)

¯ S and M ¯ i are full-symmetric in-plane base tensors defined as a function of the where M ij in-plane components of the in-plane principal directions of deformation ¯ S = 1 (n ¯2 + n ¯1 ⊗ n ¯ 2) M 12 2 ¯1 ⊗ n ¯i = n ¯i ⊗ n ¯i M Appendix B. Plane stress constraint at Gauss point level with nested loop

45

(A.4) (A.5)

1:

2: 3:

4:

5: 6:

7:

Set initial value for the logarithmic trial thickness strain tr E33 = t E33 and update the deformation component F33 = exp(λ3 ) Apply the 3D integration algorithm to obtain the stress tensor t+∆t T and the consistent tangent modulus tensor t+∆t A Check the plane stress convergence IF |t+∆t T33 | < T OL THEN Go to step 6 and EXIT. ELSE Proceed to step 4. END Update the logarithmic thickness strain by using the Newton - Raphson correction  −1 ∂2Ψ t+∆t (j) tr t+∆t (j) E33 = E33 − T33 ∂E33 ∂E33 Go to step 1. Update the in-plane consistent tangent modulus  −1 ∂2Ψ ∂2Ψ ∂2Ψ ∂2Ψ t+∆t Aα = − : : ∂Eα ∂Eα ∂Eα ∂E3 ∂E33 ∂E33 ∂E3 ∂Eα Apply kinematic mappings to compute t+∆t Sα and t+∆t Cα

46

• • • • •

Plane-stress projected large-strain multiplicative hyperelasto-plasticity The formulation includes nonlinear kinematic hardening using only Kröner-Lee decompositions Any hyperelastic function may be used for any energetic term Formulation based on elastic corrector rates, is valid for soft materials Algorithm has the same structure as a small strains 3D formulation