Fully implicit numerical integration of the Yoshida-Uemori two-surface plasticity model with isotropic hardening stagnation

The paper deals with the numerical investigation and implementation of the two-surface plasticity model (or bounding surface model). This plasticity theory allows to describe the deformation behavior under large strain cyclic plasticity and the material stress-strain responses at small-scale re-yielding after large pre-straining. A novel strategy to model the isotropic hardening stagnation is developed within a fully implicit integration scheme in order to speed up the computation and to improve the material description.


INTRODUCTION
lastic deformation on metallic materials represents an important aspect in civil engineering and in several industrial sectors (automotive, construction, naval, etc.). An underestimation of the material capability of bearing loads, or an incorrect design of components can lead to a catastrophic outcome with severe consequences in terms of lives or economic impact. This problem has been deeply investigated by many authors, resulting in a high number of analytical and numerical models that try to catch a realistic description of irreversible deformation in structural analyses and design processes. This phenomenon also occurs in sheet metal forming: the shape of a part changes during removal from the tooling due to the recovery of elastic deformation. In this case, it is important to predict accurately the springback and to compensate it correctly during the die design. In cyclic mobility problems, for example in the field of sheet metal forming, the Bauschinger P effect should be taken into account for a better prediction of springback (Uemori et al. [1])). In many typical large strain applications, the springback is a process of small-scale re-yielding after a large pre-strain. Therefore, it is necessary to define a constitutive model able to describe both the deformation behavior in large strain plasticity and the stress-strain response at small-scale re-yielding. Yoshida and Uemori [2] formulated an elastoplastic model for the description of metals behavior under cyclic loading conditions based on experimental observations on mild and high strength steels [3]. Among others constitutive theories (e.g. [4][5][6][7][8]) the Yoshida and Uemori approach appears to be the most consistent and capable to model the transient Bauschinger effect, permanent softening and the workhardening stagnation under large elastoplastic deformation. On the other hand, the increase of the computational power incentivized the research on mechanical models. The great interest around numerical simulations derives from the possibility of conducting parametric studies on components, varying geometrical features or loading conditions without additional experimental costs. Numerical codes can be applied for the simulation of service life of large structures such as bridges, skyscrapers, ships, for which full-scale experimental approaches would be unfeasible or too costly. In particular, the recurs to implicit integration schemes of material constitutive equations allows to furtherly improve the computation time, without loss of accuracy. Ghaei and Green [9] and Ghaei et al. [10] proposed a fully implicit and a semi-implicit algorithms for the integration of the two-surface model, investigating the predictive ability of the constitutive model in metal forming processes. Moreover, they enriched the formulation of the Yoshida-Uemori model by taking into account an anisotropic yield criterion, the Yld2000-2D yield function [11]. An interesting aspect in [9] and [10] is the resolution strategy for the implicit update of the workhardening stagnation, necessary to model the stress plateau observed in metals sheets subjected to cyclic loading conditions. The same numerical scheme for the workhardening stagnation was also adopted by Jia [12] to simulate the metal behavior under cyclic loading and large plastic strain. A brief description of the algorithm proposed by Ghaei and Green is presented in the following section 3.2, highlighting the main aspects and assumptions. In detail, it will be shown in section 4 that the aforementioned strategy is characterized by a non-negligible inaccuracy in the update of the workhardening stagnation surface that accumulates during loading cycles, leading to an imprecise description of the material behavior. The goal of this work is to develop a fully implicit integration scheme of the two-surface plasticity model proposed by Yoshida and Uemori, formulating a novel algorithm for the evaluation of the workhardening stagnation with an accuracy imposed by the user. The paper is organized as follows. Section 2 deals with some preliminary aspects, defining the kinematic framework accounting for large plastic deformation in metals and introducing the two-surface model. Section 3 is dedicated to describe the fully implicit integration scheme, focusing on the numerical strategy for the update of the workhardening stagnation. Both the Ghaei and Green and the proposed algorithms are presented. Section 4 reports the results obtained in the finite element analyses (FEA), showing the advantages of adopting the novel integration strategy. The limitations and future works are discussed in section 5. Lastly, conclusions are reported in section 6.

PRELIMINARIES
he two-surface model with non-isotropic hardening memory surface has been widely used to perform numerical simulation on cyclic mobility and metal forming problems [9,10,12,13]. Some limitations were pointed out by [14] in case of large plastic strains (> 20%). In the original formulation, Yoshida and Uemori expressed the constitutive model under hypoelastic-based plasticity, adopting the Jaumann spin for the definition of the co-rotational framework. This choice has been often criticized since it leads to abnormal oscillatory stress in classical shear tests (e.g. in [15], among others). The aforementioned shortcoming is discussed in the following subsection 2.1. On the other hand, the adoption of a hypobased-plasticity framework has the benefit that the constitutive equations developed for small deformation material models can be reused for large deformation simply defining a co-rotational framework. Therefore, the use of a different corotational spin does not alter the two-surface model constitutive equations, briefly recalled in subsection 2.2.

Kinematic framework
Two recent works from Jiao and Fish [16,17] investigated extensively the role of the co-rotational framework, pointing out the drawbacks of the most common objective stress rates (i.e. Jaumann, Green-Naghdi and logarithmic). In addition, Jiao and Fish proposed a new co-rotational spin, the kinetic logarithmic spin, demonstrating that this new co-rotational spin can bridge multiplicative hyper-elasto-plasticity and additive hypo-elasto-plasticity models. T The current formulation of the two-surface model adopts an additive decomposition of the strain rate D into elastic D e and plastic D p parts. The expression of the co-rotational rate of the Cauchy stress σ (plastic incompressibility is assumed  , σ σ τ Kirchhoff stress) can be written as follows: where  σ is the co-rotational kinetic logarithmic stress rate, log k Ω is the kinetic logarithmic spin, W is the continuum spin, W p is the plastic spin tensor defined by means of the material constant  according to Zbib and Aifantis [18].
is a skew-symmetric second-order tensor valued function dependent on the total strain rate D as well as on the kinetic left Cauchy-Green deformation tensor k B . This work does not consider plastic anisotropy, therefore, the contribution of the plastic spin tensor is neglected by imposing   0 . Future works will consider this aspect.

The Yoshida-Uemori model
Yoshida et al. [3] conducted a series of in-plane cyclic tension-compression tests on steel sheets in order to characterize the material deformation at large strain. Based on the observation in the experimental campaign, they formulated a constitutive model (Yoshida and Uemori [2]) to point out and to correct the shortcomings of classical models with mix isotropic and kinematic hardening. The main features of the two-surface model address three main aspects: • the reverse deformation is characterized by an initial early yielding and smooth elastoplastic transition with a rapid change of the workhardening rate followed by a softening region. • In case of mild steel sheets, the shape of the reversed stress-strain curve needs an ad hoc modeling. • The cyclic stress amplitude strongly depends on cyclic strain ranges and mean strains.
where the symbol '°' indicates the co-rotational rate,  σ is the deviatoric part of the Cauchy stress, R is the isotropic hardening function of the bounding surface, α and β are the back stress of the yield and bounding surface. B, Y, m, b, C and Rsat are material parameters,  is the plastic multiplier and H is the cumulative plastic strain. A detailed explanation of the model can be found in Yoshida and Uemori [2]. As it can be seen in Eqs. (2) 1 and (2) 2 the original model considers a classical von Mises plastic potential. Moreover, Yoshida and Uemori considered an empirical approach for the reduction of the elastic modulus E in relation to the generation of plastic deformation.
Where E0 and Ea are the Young's modulus of the virgin and infinitely large pre-strained material, respectively, and  is a material constant. The same approach has been used in this work; however, the reduction of the elastic stiffness could be better described within the continuum damage mechanics framework [19][20][21].
One of the main features of the two-surface model is the possibility of describing the stress plateau observed after the reverse loading or re-loading in the experiments. From a physical point of view, the phenomenon is related to the dissolution of dislocation cell walls during a reverse deformation. To model this behavior, Yoshida and Uemori introduced an additional surface in the stress space (see Eqn. (2)3) that can translate and expand by means of the evolution of the back stress tensor q and the radius r. In detail, the isotropic expansion of the bounding surface is allowed only if the following conditions are satisfied simultaneously: , , : : 0 The first condition in Eqn. (4)1 is fulfilled whenever the back stress of the bounding surface lies on  g ,whereas the second in Eqn. (4)2 is satisfied when the co-rotational rate of β is directed outside the non-IH surface. If both Eqs. (4)1 and (4) Therefore, the co-rotational rate of the back stress can be re-written as: ) is a material constant set by the user. Briefly, larger values of h give a rapid expansion of the non-IH surface with the consequent smaller cyclic hardening of the material. Eqn. (5) 1 is obtained directly from the consistency condition Eqn. (2)1. On the contrary, in case one or both the requirements are not satisfied, the isotropic hardening term of the bounding surface is not updated (i.e.,   0 R ). Only the back stress q is updated through Eqn. (6) where the back stress of the yield surface α can be simply obtained by means of Eqn.
(2) 8 once the terms * α and β are computed at the n+1 time step. In [2] the constitutive equations of the two-surface model are integrated by means of an explicit scheme, that requires to proceed with small time increments t to fulfill the convergence and to obtain accurate solutions. In this work the set of equations in Eqn. (2) is integrated with a fully implicit return map algorithm in the form of the closest point projection method [22,23] (CPPM). For sake of simplicity, the following subsection 3.1 shows the integration algorithm without accounting for the evolution of the non-IH surface. The workhardening stagnation requires an additional numerical scheme that is described in subsection 3.3. Lastly, subsection 3.4 summarizes the complete CPPM used for the update of all the unknowns in Eqn. (2).

CPPM without stagnation of the isotropic hardening
Firstly, the trial elastic is performed by computing the new stress state: The trial elastic stress state is used to judge if a plastic correction is required or not. In case the trial stress state is not an admissible stress state for the material (i.e.,  0 f ) a plastic correction is performed through a series of k sub-iterations It is therefore possible to write a system of five equations associated with the variables as in Eqn. (10) and to compute the solution at the k+1 sub-iteration by means of Eqn. (11).
where A is the matrix of the partial derivatives of the equations in B against the unknowns. The idea is to perform a multiequation Newton-Raphson scheme that guarantees a quadratic rate of convergence, speeding up the computation time compared to an explicit integration scheme. Moreover, the algorithm is stable and accurate for large D inputs. The subiterations k stop whenever the following condition is fulfilled: where tol1 is a threshold imposed by the user, usually 10 -8 , as in the present work.
Another benefit of the current approach consists in the possibility of computing the consistent tangent operator  TO directly after the local equilibrium is fulfilled (i.e., Eqn. (12) is satisfied) as: is a matrix with the same dimensions as elastic constant matrix E and the first term is located in the first row and second column of . A detailed explanation of the algorithm, together with iso-error maps for the evaluation of the finite step accuracy, can be found in [24,25].

Ghaei and Green algorithm for the definition of the workhardening stagnation
In large finite steps integration, the conditions in Eqn. (4) can be rewritten as the condition in Eqn. (14), which implies that the update of the workhardening stagnation has to be performed whenever the back stress  lies outside the non-IH estimated at the n step. Ghaei and Green proposed a single scalar equation for the update of  g . Assuming to know the back stress   1 1 k n β of the bounding surface, the radius of the non-IH surface at the n+1 step is approximated with a Taylor expansion as: The only unknown in Eqn. (15) is represented by the term   that can be obtained substituting the expressions in Eqn. (15) into the consistency condition of Eqn. (4)1. After mathematical manipulations (details can be found in [9])   can be computed as follows:

Novel algorithm for the definition of the workhardening stagnation
The novel approach consists in the observation that, if the condition in Eqn. (14) is fulfilled, the back stress of the bounding surface 1 n β , computed at the k+1 sub-iteration of the CPPM, should lies on the non-IH surface. Therefore, it is possible to estimate the updated  g by means of a series of additional sub-iteration j within the same k+1 sub-step. In particular, can be obtained from the previous   , 1 j n g by a Taylor expansion that neglects the quadratic terms: Figure 3: schematic representation of the evolution of the non-IH surface.
The only terms that need to be estimated are the increments of the back stress   are updated through Eqn. (19). The algorithm stops whenever the following conditions is fulfilled: From a theoretical point of view, the scheme for the update of the non-IH surface can be seen as a sort of return mapping where, instead of correcting the back stress   1 1 k n β , the values of the radius and the back stress of  g are updated (see Fig. 3). It is important to highlight that the condition in Eqn. (21) expresses exactly the updated non-IH surface that passes through the current value of   1 1 k n β . Moreover, the threshold tol2 can be set independently form tol1, depending on the accuracy required. In this work both tol1 and tol2 were both set equal to 10 -8 .

CPPM with stagnation of the isotropic hardening
This section summarizes the steps necessary for the computation of all the variables in Eqn. (7) from the step n to the subsequent step n+1. The procedure is reported in the flow chart of Fig. 4. Subsequently, the isotropic hardening term R for the bounding surface is updated. Instead, in case the condition in Eqn. (14) is false, the non-IH surface is not updated and is updated,   1 n n r r . 5. Lastly, the fulfillment of the local convergence at the n+1 step and k+1 subiterations is checked with Eqn. (12). If Eqn. (12) is true the algorithm proceeds to the next time step, otherwise it is necessary to repeat the procedure from point 2.

NUMERICAL ANALYSES
he constitutive equations presented in the previous subsections 3.1and 3.3 were implemented via user subroutine for the commercial code Abaqus (ver. 6.14-4). The present section deals with the results of the numerical analyses divided into two parts. The first part aims to show the benefit of the novel approach for the update of the non-IH surface compared with the Ghaei and Green algorithm. The second part shows the ability of the model to reproduce the experimental results carried out by Yoshida et al. [3] on a mild (i.e., SPCC) and high strength (i.e., SPFC) steel sheets. The material parameters are reported in Tab. 1, while the geometry of the sample and the sketch of the mesh and boundary conditions are reported in Fig. 5a and b, respectively. The FEA were carried out considering 1695 hexahedral elements with reduced integration (i.e., Abaqus element C3D8R) and applying fully reversed displacement conditions on the top of the sample. To reduce the computation time only 1/8 of the sample was considered (see red dashed line in Fig. 5a). Moreover, to reproduce the experimental set up, the loading condition was controlled by a sensor positioned to simulate the strain gauge on the specimen. Whenever the nominal strain on the sensor exceeded the ±5% axial strain the direction of the load was inverted (e.g., from tensile to compressive and vice versa).

Ghaei and Green algorithm vs novel approach
The purpose of this section is to compare the Ghaei and Green formulation against the proposed approach. The sample in Fig. 5b is subjected to three fully reversed loading cycles; the material parameters correspond to the SPCC steel and the parameter h is set to 0.9 to enhance the expansion of the radius r over the back stress q. The results by the explicit integration in Yoshida and Uemori [2] are indicated with black solid lines, the FEA carried out with the Ghaei and Green algorithm with blue markers while the current approach is displayed with red markers.
As it can be seen in Fig. 6a, the proposed approach gives a better description of the material behavior catching a realistic stress plateau during the first reverse loading. The Ghaei and Green formulation overestimates the stress plateau resulting in a less accurate description of the stress-strain curves in the remaining cycles. The performance of the two formulations can be better observed in the graphs of Fig. 6b and c, reporting the evolution of q and r against the nominal strain. While the current approach is able to describe the evolution of the back stress and the radius with good approximation. The algorithm by Ghaei and Green suffers from an underestimation of the back stress evolution and an initial overestimation of the radius. This results can be explained by the fact that r in Eqn. (15)

Numerical results on SPCC and SPFC steel sheets
The aim of this section is to show the performance of the fully implicit integration scheme against the experimental results carried out by Yoshida et al. [2]. Fig. 7a and b report with black solid lines the stress-strain responses during cyclic loading with total strain ranges of 4% and 10% for the mild and high strength steel, respectively. The red lines show the stress-strain curves obtained with the numerical model. As it can be seen the FE simulations can catch quite well the material response, especially the workhardening stagnation during the reverse loading, proving the correct implementation of the proposed approach. In general, the stress plateau is observed in the mild steel, whereas the high strength steel gives a hardening response. The elastic stiffness seems slightly underestimated in the numerical analyses, probably due to the use of the sensor for the control of the loading conditions. Overall, the numerical results are in good agreement with the experimental data.

DISCUSSION
he benefit of the Ghaei and Green formulation lies on the use of a single scalar equation for the update of the non-IH surface. On the other hand, the proposed algorithm requires an iterative procedure for modeling of the workhardening stagnation. However, it has been shown how the Ghaei and Green algorithm is characterized by less accuracy due to the lack of control on the mutual position between the final non-IH surface and the bounding surface back stress. The novel approach can guarantee the fulfillment of the consistency condition for the non-IH with an accuracy set by the user. Moreover, Fig. 7c and d report the convergence graphs of the CPPM and of the update algorithm for the workhardening stagnation, respectively. The data are reported for the quadrature point of the element experiencing the largest amount of cumulative plastic strain H (the curves are reported at 10-50-100% cumulative plastic strain values). As it can be seen the T novel approach does not require a large amount of sub-iteration j (Fig. 7d) and it does not influence the convergence rate of the k sub-iterations (Fig. 7c).

CONCLUSIONS
he work presented a novel algorithm for the computation of the evolution of the workhardening stagnation surface in the Yoshida-Uemori model. The numerical scheme guarantees a quadratic rate of convergence for the local and global equilibrium without loss of accuracy even at large integration steps. The current implementation does not consider an anisotropic yield criterion, but it assumes a classical von Mises plastic potential. In case of hot-rolled or cold-rolled metals sheets, anisotropy plays a fundamental role in the material behavior and therefore a J2 plasticity might lead to inaccurate results. Future works will consider this aspect. A recent work from Xie et al. [26] formulated the workhardening stagnation phenomenon with a non-IH defined in the strain space rather than in the stress space, with an additional recovery term. Xie's model showed a better description of the hysteresis loops under various amplitudes compared with the original two-surface model. Future works will investigate this aspect. It is worth mentioning that a definition of the workhardening stagnation in the strain space does not invalidate the proposed approach, that can be still adopted. Lastly, it should be pointed out that the level of deformation in the current work is quite limited. The results carried out with the Jaumann co-rotational rate and the kinetic logarithmic spin showed negligible differences. Future works will consider higher deformation regimes.