State update algorithm for associative elastic-plastic pressure-insensitive materials by incremental energy minimization

This work presents a new state update algorithm for small-strain associative elastic-plastic constitutive models, treating in a unified manner a wide class of deviatoric yield functions with linear or nonlinear strain-hardening. The algorithm is based on an incremental energy minimization approach, in the framework of generalized standard materials with convex free energy and dissipation potential. An efficient method for the computation of the latter, its gradient and its Hessian is provided, using Haigh-Westergaard stress invariants. Numerical results on a single material point loading history and finite element simulations are reported to prove the effectiveness and the versatility of the method. Its merit turns out to be complementary to the classical return map strategy, because no convergence difficulties arise if the stress is close to high curvature points of the yield surface.


INTRODUCTION
he solution of inelastic structural boundary value problems in a typical finite element implementation requires the integration of the material constitutive law at each Gauss point.The complexity of hardening elastic-plastic constitutive equations requires their time discretization and numerical integration on a sequence of finite time steps.This procedure is referred to as the state update algorithm.The most traditional way to formulate the associative plasticity constitutive law is based on the assumption of a convex yield function and the adoption of a normality law to determine the plastic flow.A widely used class of constitutive updates exploiting this framework is the so-called return mapping algorithms.While the basic idea dates back to the work of Wilkins [1], all the return mapping algorithms share the use of an elastic-predictor, plastic-corrector strategy.In particular, in a first step the algorithm checks whether the elastic trial state is plastically admissible.If such is not the case, a solution for the time-discrete plastic evolution equations is sought for in the plastic corrector step.From a geometric standpoint, the return map algorithm amounts to finding the closest point projection, in a suitable norm, of the elastic trial state on the yield domain.To this purpose, a typical method pursued among others in [2][3][4][5][6] consists in the iterative T solution by a Newton-Raphson scheme.That strategy appears to be attractive because of the quadratic rate of convergence of Newton method, preserved by the adoption of an algorithmically consistent tangent operator.Details on the subject can be found in several recent textbooks, e.g.[7][8][9][10][11].However, in many cases of practical interest, the localconvergence characteristic of Newton schemes does not guarantee global convergence [12][13][14].For instance, convergence may not be attained for elastic trial stress in the neighborhood of high curvature points of the yield surface [15].In the literature, various techniques have been proposed to avoid this drawback, as the use of line search methods, sub-stepping, adaptive sub-stepping, or other ad-hoc integration algorithms, see for example [16][17][18][19][20][21].An alternative approach is to set the elastic-plastic problem in the framework of normal-dissipative standard media in the sense of Halphen and Nguyen [22].The evolution of internal variables is assumed to be governed by a convex scalar dissipation function which is the support function of the elastic domain (e.g., [23]).Consequently, the flow law results in a statement of the classical postulate of maximum plastic dissipation [7].Exploiting this formulation, several state update algorithms have been proposed, either at local level [24][25][26] or at finite element level [27,28].Conceptually, the underlying basic idea is that the evolution of internal variables in a finite time step incrementally minimizes a suitable convex functional, given by the sum of the internal energy and the dissipation function.This paper focuses on the incremental energy formulation, for elastic-plastic hardening materials characterized by isotropic deviatoric yield function and associative flow law, in the framework of infinitesimal deformation.A two-step algorithm is proposed to perform the material state update.In the first step, an elastic prediction of the updated material state is carried out.In case it is not plastically admissible, the Newton-Raphson method is adopted to solve the constitutive variational problem.An efficient strategy to compute the dissipation function is proposed.In particular, adopting the Haigh-Westergaard representation (e.g.see [29][30][31]), the problem is reduced to a non-linear scalar equation.Moreover closed-form expressions of the gradient and of the Hessian of Haigh-Westergaard coordinates, as well as of the dissipation function, are presented.With the aim of proving the robustness and stability of the proposed algorithm, numerical tests on a single integration point and FE simulations are provided.This numerical approach appears to be complementary to the classical return map strategy, because no convergence difficulties arise if the stress is close to points of the yield surface with high curvature.The paper is organized as follows: a first part is concerned with a brief discussion on the incremental energy minimization in hardening elastoplasticity.Then the computation algorithm of dissipation function and the proposed state update algorithm are extensively treated.Finally a selected number of numerical tests are reported.
, in practice a Gauss point in a typical finite element solution.In the case of infinitesimal deformation, the total strain ε , i.e. the symmetric part of the displacement gradient, is additively decomposed as: where e ε and p ε are the elastic and plastic parts, respectively.The kinematic description of the material state is completed by a set of strain-like internal variables  ε α Assuming the latter to be strictly convex, the constitutive equations follow , , , where (•)  is the sub-differential operator of non-smooth convex functions (see [22,32] and references therein).In the context of elastic-plastic materials, the generalized stresses σ , q are constrained to belong to an admissible set, denoted as elastic domain and introduced by means of the yield function   Hereafter we suppose the yield function to be convex; consequently the elastic domain is convex itself.

L
The model is supplemented by the evolution law of plastic strain p ε and strain-like internal variables α .With the assumption of normal-dissipative or standard material behavior, i.e. of associative plastic flow, it is customary to introduce a scalar dissipation function defined as the support function of the elastic domain and to express the evolution law in the form: It is worth noting that the dissipation function is convex, non-negative, zero only at the origin and positively homogeneous of degree one; moreover it is nondifferentiable at the origin, where its sub-differential set coincides with the elastic domain [7].
Substituting the decomposition (1) in the constitutive law (2) and adding the result to (5), the constitutive differential equations are obtained: , , often referred to as Biot's equation of standard dissipative systems, see [22,24,33] for instance.We now proceed with the numerical approximation of the rate Eq. ( 6).In particular, we adopt a backward Euler integration scheme and assume that plastic strain p ε and strain-like internal variables α vary linearly in each time step 1 ] [ , n n t t  .Consequently, by exploiting the degree-one positive homogeneity of the dissipation function, the incremental form of (6) follows: with the notation     Eq. ( 7) represent the Euler-Lagrange first order conditions associated to the incremental minimization problem where e,trial α define the elastic trial state, i.e. the material state corresponding to no plastic evolution.Therefore, the incremental minimization problem (8) determines the internal state of the material for finite increments of time.We remark that the present formulation is similar to the one proposed in [24], under the assumption of rateindependent evolution and linear variation in time of plastic strain p ε and strain-like internal variables α .

COMPUTATION OF THE DISSIPATION FUNCTION
inematic and isotropic hardening are distinguished by assuming that the yield function can be represented as , where k q [resp., i q ] is the stress-like kinematical [resp., isotropic] hardening variable.Accordingly, the dissipation function is given by where k α [resp., i   ] is the strain-like kinematical [resp., isotropic] hardening variable.Setting With the assumption that the yield function is isotropic, it can be represented, with a slight abuse of notation, as are the Haigh-Westergaard coordinates of σ (see Appendix A).Since the maximum of the scalar product of two symmetric tensors is attained when they are coaxial (see, e.g., [34]), in (12) it is possible to assume that σ is coaxial with p ε and to express their scalar product making use of the Haigh-Westergaard representation.Accordingly, the dissipation function is given by

Deviatoric yield function, perfect plasticity or kinematic hardening
The yield function is assumed to be of the form where g is a 2 C , positive, even, periodic function with period 2 3  can be interpreted as the initial yield limit in tension if   According to (14), the dissipation function is given by This implies: where or, equivalently, . Because the supremum is attained on the boundary of that set, a simple computation yields which can be recast as: (20) This is a scalar equation in the unknown . By applying Dini's theorem, its first derivative turns out to be: (21) where g and its derivatives are computed at   by (20) it follows that whence, recalling ( 18) and ( 21 where the derivatives on the right-hand sides of ( 23 Simple expressions for the gradient and the Hessian of p  ε and p   are given in Appendix A. The same argument holds if kinematic hardening is present.

Deviatoric yield function, isotropic hardening
The yield function is assumed to be whence the dissipation function follows: This implies: Accordingly, the dissipation function turns out to be: whence it is finite provided that In case of hardening material, recalling (25), the solution 0 i y q   of the supremum in ( 27) is unfeasible.Hence it can be assumed that ( 16) prevails, complemented with The same treatment holds if also kinematic hardening is present.

STATE UPDATE ALGORITHM
n this section, under the assumption of deviatoric yield function, we introduce a reduced formulation for the incremental minimization problem (8) and discuss the state update algorithm proposed for its solution.The following common choice for the free energy function e ( , ) where W and  , respectively, denote the elastic and hardening potentials, and may account for nonlinear elastic and/or hardening behavior.The contributions from kinematic and isotropic hardening are distinguished in the latter, as follows: As a consequence, substituting k α [resp., i   ] from (11) [resp.(29)], the incremental energy problem ( 8) is transformed into: Under isotropy assumption, the elastic potential is decomposed as: where sph W and dev W denote the spherical and deviatoric parts of W , e e denotes the deviatoric part of e ε , and tr denotes the trace operator.Moreover, without loss of generality, the kinematic hardening potential is assumed to depend only on the deviatoric part of k .
α Therefore the incremental energy problem (32) is further reduced to: where p e is the deviatoric part of p ε .Recalling that the dissipation function   p D e is singular at the origin p   e 0, the state update algorithm first computes an elastic predictor.In case the corresponding state is plastically admissible, the step is elastic and the infimum in ( 34) is attained at the origin.Conversely, having excluded the singularity, a smooth incremental energy minimization problem can be solved.By introducing the linear projector operator on the deviatoric tensor subspace dev / 3    I I   , the Euler-Lagrange equation of ( 34) is: where the gradients of dev W , k  , and D are taken, respectively, with respect to e e , k α , and p e .Eq. ( 35) is solved by the Newton-Raphson method, using the consistent tangent:   35), the algorithmic consistent elastoplastic tangent is given by: Remark.The present formulation, when specialized to particular choices of yield function and elastic/hardening potentials, gets simplified.As an example, we consider a linear elastic and hardening behavior with von Mises criterion.The Helmoltz free energy and the hardening potentials are given by: the yield function and the relevant dissipation function are given by: The Euler-Lagrange Eq. ( 35) becomes: Introducing the deviatoric trial back-stress trial e,trial trial noting that the latter results to be parallel to ) As a consequence, the well-known expression (see, e.g., [35]) for the algorithmic consistent elastoplastic tangent results: NUMERICAL TESTS n this section we assess the numerical reliability and robustness of the proposed state update algorithm.Its convergence properties are explored by a single increment test and compared with the performances of the classical return mapping approach.The numerical results prove that the two strategies have complementary convergence properties.Then, we consider the integration of a non-proportional stress-strain loading history of a material point, and finite element simulations.

Single increment test
In [17] Armero and Pérez-Foguet assess the converge properties of return mapping algorithms by using the following single increment test.Assuming the initial state to be virgin, a total strain ε (coinciding with the elastic trial strain e,trial ε ) is applied to the Gauss point.The total strain is defined in terms of its Haigh-Westergaard coordinates in the deviatoric

I
Fig. 1 compares the performances of the return mapping and incremental energy minimization algorithms in terms of number of iterations needed for convergence.In Fig. 1(a) it is clear that when the elastic trial state in stress space is close to high curvature points of the yield surface, i.e. for 0 , the return mapping algorithm has convergence difficulty.As depicted in Fig. 1(b), for the same elastic trial state the proposed state update algorithm needs 1-3 iterations to reach convergence.A converse situation happens around 30 . This result is explained observing that high curvature points in the yield surface correspond to low curvature points in the graph of the dissipation function, and viceversa (see Appendix B).

Pointwise mixed stress-strain test
In order to explore the algorithm effectiveness, we consider a pointwise mixed stress-strain test as reported in [36].In particular, a non-proportional stress-strain history is adopted.respectively.In Fig. 2 the numerical solution, expressed in terms of the stress-strain history, proves to be in perfect agreement with the analytical one (for example, see [35]).

Stretching of a perforated rectangular plate
A classical benchmark for the implementation of plasticity models is the stretching of a thin perforated plate along its longituindal axis (see [10], for instance).This problem is solved under the usual assumption of plane stress, which is straightforwardly imposed within the present formulation.To this end, it suffices to minimize the incremental energy functional in (34) also with respect to the total strain component normal to the stress plane.The plate is 20 mm wide, 36 mm high, 1 mm thick and presents a circular hole of radius 10 mm at its center.The geometry and boundary conditions are depicted in Fig. 3(a).The material is supposed to be isotropic linear elastic, with Young modulus 70 GPa E  and Poisson coefficient 0.2   .A von Mises yield criterion (see Appendix B), with yield stress y0 0.243 GPa   and linear isotropic hardening of modulus i 0.2 GPa H  , is adopted.
Because of the symmetry, only one quarter of the plate is considered in the analysis by applying the appropriate boundary condition on the symmetry edges.The finite element mesh adopted is constituted by 576 OPT triangular elements (see [38], for instance), with a total number of 325 nodes.The analysis is conducted under displacement control, assigning a vertical displacement on the top constrained edge.In Fig. 3(b) the diagram of the total reaction force versus the prescribed displacement is reported.The numerical solution is in good agreement with the reference one [10].In Fig. 3(c [39]).

Pinching of a cylindrical shell with diaphragms
In this section we consider the application of the proposed state update algorithm to a geometrically nonlinear problem.
In particular, we refer to the problem of a cylindrical shell, constrained by rigid diaphragms at end-sections, and pinched by concentrated forces in the mid-section [39].The cylinder has radius of 300 consistent units (c.u.), while its total length is of 600 c.u.The initial shell thickness is 3 c.u. Fig. 4(a) shows the geometry and boundary conditions (we assume that the presence of rigid diaphragms constrains the in-plane degrees of freedom of the cylinder end sections).We assume the material to be isotropic linear elastic, with Young modulus 3000 E  c.u. and Poisson coefficient 0.3   .A von Mises yield criterion (see Appendix B), with yield stress y0 24.3   c.u. and linear isotropic hardening of modulus i 300 H  c.u., is adopted, under plane-stress assumption.Due to symmetry, only one-eighth of the cylinder is considered in the analysis by applying appropriate boundary conditions on the symmetry edges.The geometrical nonlinearity is treated by a corotational approach [40,41], and the finite element mesh is constituted by 24×24 OPT-DKT triangular elements (see [38,42], for instance).The analysis is conducted under displacement control, up to a total displacement of 300 c.u., coinciding with the cylinder radius.Fig. 4(b) depicts the deformed configuration at the final load level, while in Fig. 4(c) load-displacement diagram is plotted.
The numerical solution appears to be in good agreement with the reference one [39].The slight difference at high displacement values may be attributed to different finite-element types used in the computations.

CONCLUSIONS
new state update algorithm for small-strain associative elastic-plastic constitutive models was presented.It is set within the theoretical framework of generalized standard materials with internal variables, under the assumption of rate independence.The evolution of internal variables in a finite time step obeys an incremental minimization principle involving a suitable convex functional, given by the sum of internal energy and dissipation function.An efficient strategy to compute the latter and its gradient and Hessian is proposed, using Haigh-Westergaard stress invariants, for an ample class of isotropic deviatoric yield functions with linear or nonlinear strain-hardening.Numerical results proved the effectiveness and the versatility of the methodology on a single material point state update, for a given loading history under mixed stress/strain control, and showed its reliability as an efficient constitutive solver in conjunction with finite element simulations.The proposed algorithm turned out to be complementary to the classical return mapping strategy, because no convergence difficulties arise if the stress is close to points of the yield surface with high curvature.Hence it may improve the robustness of available plastic solvers through a hybrid approach.Future work will involve the modeling of pressure-sensitive yield surfaces and non-associative plastic flow rules, and the testing of the solution algorithm on extensive large-scale FE simulations of engineering-relevant problems.
where ε is a symmetric second-order tensor and e is its deviatoric part , defined by: I is the second-order identity tensor,  is the usual dyadic product between second-order tensors, and  is the fourthorder identity on second-order symmetric tensors.The gradients and the Hessians of those invariants are given by: Here (•)  denotes the derivative with respect to ε ,  is the fourth-order null tensor,  denote the square tensor product between second-order tensors [43].The analysis presented herein is however greatly simplified by using another set of invariants, known as Haigh-Westergaard coordinates [29,44], defined by: The gradients and Hessians of  and  are straightforwardly derived: whereas the derivation of   and   is somewhat more involved.To this end, the argument in [45] is briefly sketched, and the trihedron 1 2 3 { , , } f f f of second-order symmetric tensors is introduced: Using (A4), (A5) and (A3), it turns out that: The latter equation yields by (A3) and (A4): whence   follows by (A6).On the other hand, Eq. (A7) can be recast as and, by differentiating with respect to ε , it leads to: Noting that by (A9), (A3) and (A7) it turns out that: and by (A6), (A5) and (A3) it turns out that: it is a simple matter to verify that (A10) yields: Alternative expressions for   and   , more effective from a computational point of view, are: The expressions for   and   derived in (A8), (A6) and (A13), or the alternative expressions in (A15), get greatly simplified when computed in an orthonormal basis of eigenvectors of ε .The tensor 1 f is spherical, thus it is represented by   In passing, it is noted that  may be assumed to belong to [0, 3]  , implying that the eigenvalues of e are sorted in descending order.Substituting 1 f and 2 f into (A8), it turns out that 3 f is represented by: Analogously, substituting the above representations into (A13), the representation of   is obtained.The only nonzero terms turn out to be: Expressions (A18) was reported in [45]; to the best of Authors' knowledge, (A19) is new.It is emphasized that the singularity of   at integer multiples of 3  , due to the sin3 term in the denominator of (A8) or (A15), turns out to be removable, because it disappears in (A18).On the other hand, the "shear components" of   are indeed singular at integer multiples of 3  , as apparent by (A19).However,   only appears in (24), multiplied by D  .Using (23) and recalling that ' g vanishes at integer multiple of / 3  , it turns out that the product D    has removable singularities there, thus allowing for a stable numerical computation of D  .

APPENDIX B: GALLERY OF DEVIATORIC YIELD FUNCTIONS
Von Mises yield function eferring to the form reported in ( 14) and (25), the von Mises yield function is obtained by assuming (see [45], for instance):

Von Mises-Tresca yield function
The von Mises-Tresca yield function corresponds to the choice (see [17] and references therein): The corresponding conjugated stress-like variables are the stress n   σ  and the stress-like internal variables n   q  , respectively.By standard thermodynamic arguments, the stress-like variables are derived from an Helmholtz free energy function   e , .


and its first and second derivatives with respect to p   easily follow: ) are performed with respect to  σ and computed at   p    σ ε .Finally the gradient and the Hessian of   p D ε are straightforwardly obtained from (16): (34).Exploiting the degree-one positive homogeneity of the dissipation function, the Taylor expansion about condition results in a scalar algebraic equation for the Lode angle p

Figure 2 :Figure 3 :
Figure 2: Pointwise mixed stress-strain test.Comparison between analytic and numerical solution for von Mises yield surface.

Figure 3 :
Figure 3: Pinching of a cylindrical shell with diaphragms.(a) Geometry, boundary conditions and finite element mesh; (b) Deformed configuration; (c) Comparison of load-displacement diagram between the proposed algorithm and reference solution (see[39]).

2 f
the diagonal matrix with the enclosed eigenvalues.Recalling (A7), 2  where i , 1,..., 3 i e  are the eigenvalues of e , which are the roots of the characteristic polynomial of e : ), it is a simple matter to verify that those roots are is represented by:

Fig. 5
Fig.5shows the corresponding yield surface and dissipation function graph.It is worth noting that the dissipation function results convex, non-negative, zero only at the origin and positively homogeneous of degree one.The latter feature is clear from its cone-like graph.
m determines the shape of the yield surface.In particular in the limit of 0 m  [ resp., m   ] the von Mises circle [resp.the Tresca hexagon] is recovered.High curvature of the yield surface corresponds to high values of the parameter m .For 20 m  the discrepancy between this yield surface and the Tresca hexagon is less than 2 per cent [17]: Fig. 6 depicts the von Mises-Tresca yield surface and the corresponding dissipation function graph.Fig. 6(b,d) show that high curvature points in the yield surface correspond to low curvature points in the graph of the dissipation function, and vice-versa.

Table 1 :
Pointwise mixed stress-strain test.Controlled strain history: in between two next time units, strain components vary linearly.