Microsoft Word - numero_57_art_10_3122 R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 114 Fully implicit numerical integration of the Yoshida-Uemori two-surface plasticity model with isotropic hardening stagnation Riccardo Fincato, Seiichiro Tsutsumi University of Osaka, JWRI, Japan fincato@jwri.osaka-u.ac.jp, http://orcid.org/0000-0003-4058-0460 tsutsumi@jwri.osaka-u.ac.jp Alex Zilio, Gianluca Mazzucco, Valentina Salomoni University of Padova, DICEA, Italy alex.zilio@studenti.unipd.it, gianluca.mazzucco@dicea.unipd.it, valentina.salomoni@dicea.unipd.it ABSTRACT. 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. KEYWORDS. Bounding surface model; Return mapping; Workhardening stagnation; Finite strain; Finite element analyses. Citation: Fincato, R., Tsutsumi, S., Zilio, A., Mazzucco, G., Salomoni, V., Implicit numerical integration of the Yoshida-Uemori two-surface plasticity model with isotropic hardening stagnation, Frattura ed Integrità Strutturale, 57 (2021) 114-126. Received: 19.05.2021 Accepted: 03.06.2021 Published: 01.07.2021 Copyright: © 2021 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 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 https://youtu.be/aJprgY6tBy4 R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 115 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–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 co- rotational 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 R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 116 The current formulation of the two-surface model adopts an additive decomposition of the strain rate D into elastic De and plastic Dp parts. The expression of the co-rotational rate of the Cauchy stress σ (plastic incompressibility is assumed ,σ σ τ Kirchhoff stress) can be written as follows:                   log log log log , k k k p k k p p p σ σ Ω σ σΩ Ω W W W B σ D W σD D σ (1) where σ is the co-rotational kinetic logarithmic stress rate, logkΩ is the kinetic logarithmic spin, W is the continuum spin, Wp is the plastic spin tensor defined by means of the material constant  according to Zbib and Aifantis [18]. logkW 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 kB . 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. The model was formulated accounting for two separate surfaces: the yield surface f (Eqn. (2)1), that considers only kinematic hardening (Eqn. (2)8), and the bounding surface F (Eqn. (2)2), characterized by both isotropic and kinematic hardening (Eqs. (2)9 and (2)7). Additionally, the workhardening stagnation can be described by adopting a third surface, the non-isotropic hardening surface g (non-IH hereafter) (Eqn. (2)3). A novel algorithm for the numerical definition of the non-IH surface is the object of this work. Briefly, the constitutive equation of the two-surface model are presented in Eqn. (2), referring to Fig. 1. Figure 1: Sketch of the two-surface model.                                              * * * * * * * * * ( ) 3 2 ; 3 2 ( ) 3 2 2 3 ; 3 2 ; ; 2 3 ; ; 1 ; 2 3mHsat f Y F B R g r Ca a a B R Y f f m b R R e H dt σ α σ β β q α N N α σ α N N σ α β N β α α β (2) 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 R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 117 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.         0 0 1 H aE E E E e (3) 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–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:                    2 2, , 3 2 0 , , : : 0 g r r g r β q β q β q β β q β β (4) 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)2 are true, the bounding surface expands (i.e.,  0R ) following the evolution law in Eqn. (2)9 and the back stress q and the radius r are updated according to the following formulas: a) b) Figure 2: Sketch of the non-IH surface a) in case of non isotropic hardening (  0R ), b) in case of isotropic hardening (  0R ).                    2 3 : , 2 3 : 2 r rr r h r β q β q β q β q β (5) Therefore, the co-rotational rate of the back stress can be re-written as: R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 118        3 :1 2 h r r β q β q (6) where h (  0 1h ) 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.,  0R ). Only the back stress q is updated through Eqn. (6), while  0r . NUMERICAL INTEGRATION onsidering an equilibrium state at the generic time   0,nt t the aim of the integration algorithm is the computation of the unknowns at the subsequent time    1n nt t t with   1 0,nt t , given the total strain increment D . In detail, the set of the unknowns can be expressed as:   *, , , , , ,H rx σ α β q (7) 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:    1 : trial n nσ σ E D (8) 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.,  0f ) a plastic correction is performed through a series of k sub-iterations (starting by k = 0 and assigning the internal variable, except the stress state, as       0 1 k n n ). If the evolution of the non-IH surface is not considered, the set of the unknows in Eqn. (2) can be reduced to the following set of five unknowns:   *, , , , Hx σ α β (9) 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). C R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 119                                                             1 1 1 1 1 *, 1 1 *, 1 *, 1 *, 1 1 1 1 1 3 2 3 2 ( ) 2 3 2 3 2 3 k k n n k k k n n n k nk k k k n n n n n k k k n n n n k n n f Y F B R Ca a m b H H σ α σ β B α α N N β β N β (10)               11 1 1 1 1 1 1 1 k k k k k k n n n n n nx x A B x x (11) where A is the matrix of the partial derivatives of the equations in B against the unknowns. The idea is to perform a multi- equation 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 sub- iterations k stop whenever the following condition is fulfilled:    1 1 1 k n tolB (12) 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:      1( 1) 1 12 :TO knA E (13) where     11 1 12 k nA 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     11 1 k nA . 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  1 1 k nβ 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:           2 21 1 1 1, , 3 2 0 k k n n n n n ng r rβ q β q (14)                                              2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 3 : 1 1 k n n n n n n k k n n n n nk n n n n r r r r h χ β β β β β qχ χ β q (15) R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 120 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:                   21 1 2 1 1 2 3 3 : 3 : 4 : 2 1 2 k k n n n n n n n n h h r r χ β χ β χ χ (16) Lastly, the back stress 1nq and the radius rn+1 are updated by substituting the computed  into Eqn. (15)3 and then 1nχ in Eqn. (15)1. 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 1nβ , 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,    1 , 1 j ng can be obtained from the previous  , 1 j ng by a Taylor expansion that neglects the quadratic terms: Figure 3: schematic representation of the evolution of the non-IH surface.                                  1 1 1 , 1 , 1 1 1 11 , 1 , 1 1 with for 0 j j j j j j n n n n nn j n n j n n j n n g g g g r r g g j r r q q q q (17) The only terms that need to be estimated are the increments of the back stress  1 1 j nq and of the radius   1 1 j nr . To compute those terms the following variable  is introduced:                 1 1 1 1 1 11 1 2, 1 1 3 : 2 k j k n n nj n j nr β q β (18) From Eqs. (5) and (6) it is possible to write:                          1 1 1 1 1 1 1 1 1 1 1 1 1j j k jn n n n j j j n n n h r h r q β q (19) R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 121 At this point the only unknown of the problem is represented by   1 1 j n that can be obtained by substituting Eqn. (19)1 and 2 into the linear expansion of Eqn. (17). After manipulations it is possible to write:                         , 11 1 1 1 1 1 11 1 : j nj n j j k j j n n n nn g g g h h r r β q q (20) Once computed   1 1 j n , function of all the variables at the j sub-iteration, the values of   1 1 j nq and   1 1 j nr are updated through Eqn. (19). The algorithm stops whenever the following conditions is fulfilled:     1 , 1 2 j ng tol (21) 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. Figure 4: Flow chart of the CPPM with isotropic hardening stagnation. 1. Perform the trial elastic. 2. In case the plastic correction of the stress is needed, the scheme applies the procedure explained in subsection 3.1 through a series of k sub-iterations. 3. Within the same k sub-iteration and after the computation of the unknowns in Eqn. (11), the fulfillment of the condition in Eqn. (14) is verified. 4. The condition in Eqn. (14) is true, the algorithm updates the non-IH surface with a series of j sub-iterations by means of the procedure explained in sub-section 3.3 until Eqn. (21) is verified. 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  0R . Only 1nq is updated,  1n nr r . 5. Lastly, the fulfillment of the local convergence at the n+1 step and k+1 sub- iterations 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. R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 122 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). Material SPCC SPFC E0 [GPa] 206 206 Ea [GPa] 152 159 υ 0.3 0.3 ξ 30.8 61 Y [MPa] 124 360 Rsat [MPa] 190 143 C 500 300 b [MPa] 9 66 m 12 26 B [MPa] 168 478 Table 1: Material constants for the SPCC and SPFC steels. a) b) Figure 5: a) Geometry of the sample, the red dashed line indicates the portion modeled in the FEA (dimensions in mm). b) boundary conditions, mesh and sensor location for the control of the loading conditions. T R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 123 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)1 is evaluated with a Taylor expansion that neglects the contributions of  2r . Moreover, the condition in Eqn. (4)1 is not checked after performing the update of q and r. On the contrary, the novel approach computes simultaneously the back stress and the radius through an iterative procedure that guarantees that the final non-IH surface passes through the current back stress of the bounding surface. a) b) c) Figure 6: a) Nominal stress vs nominal strain curves. Back stress q and radius r evolution carried out with b) proposed approach, c) Ghaei and Green formulation. R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 124 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. a) b) c) d) Figure 7: Numerical vs experimental results on a) SPCC steel and b) SPFC steel. c) Local convergence of the CPPM, d) Convergence of the proposed algorithm for the update of the non-IH surface. 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 R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 125 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. REFERENCES [1] Uemori, T., Okada, T., Yoshida, F. (2000). FE Analysis of Springback in Hat-Bending with Consideration of Initial Anisotropy and the Bauschinger Effect, Key Eng. Mater., 177–180, pp. 497–502, DOI: 10.4028/www.scientific.net/KEM.177-180.497. [2] Yoshida, F., Uemori, T. (2002). A model of large-strain cyclic plasticity describing the Bauschinger effect and workhardening stagnation, Int. J. Plast., 18(5–6), pp. 661–686, DOI: 10.1016/S0749-6419(01)00050-X. [3] Yoshida, F., Uemori, T., Fujiwara, K. (2002). Elastic–plastic behavior of steel sheets under in-plane cyclic tension– compression at large strain, Int. J. Plast., 18(5–6), pp. 633–659, DOI: 10.1016/S0749-6419(01)00049-3. [4] Sanchez, L.R. (2010). Modeling of springback, strain rate and Bauschinger effects for two-dimensional steady state cyclic flow of sheet metal subjected to bending under tension, Int. J. Mech. Sci., 52(3), pp. 429–439, DOI: 10.1016/j.ijmecsci.2009.11.002. [5] Li, K.P., Carden, W.P., Wagoner, R.H. (2002). Simulation of springback, Int. J. Mech. Sci., 44(1), pp. 103–122, DOI: 10.1016/S0020-7403(01)00083-2. [6] Chun, B.K., Jinn, J.T., Lee, J.K. (2002). Modeling the Bauschinger effect for sheet metals, part I: theory, Int. J. Plast., 18(5–6), pp. 571–595, DOI: 10.1016/S0749-6419(01)00046-8. [7] Chun, B.K., Kim, H.Y., Lee, J.K. (2002). Modeling the Bauschinger effect for sheet metals, part II: applications, Int. J. Plast., 18(5–6), pp. 597–616, DOI: 10.1016/S0749-6419(01)00047-X. [8] Gau, J.-T., Kinzel, G.L. (2001). A new model for springback prediction in which the Bauschinger effect is considered, Int. J. Mech. Sci., 43(8), pp. 1813–1832, DOI: 10.1016/S0020-7403(01)00012-1. [9] Ghaei, A., Green, D.E. (2010). Numerical implementation of Yoshida–Uemori two-surface plasticity model using a fully implicit integration scheme, Comput. Mater. Sci., 48(1), pp. 195–205, DOI: 10.1016/j.commatsci.2009.12.028. [10] Ghaei, A., Green, D.E., Taherizadeh, A. (2010). Semi-implicit numerical integration of Yoshida–Uemori two-surface plasticity model, Int. J. Mech. Sci., 52(4), pp. 531–540, DOI: 10.1016/j.ijmecsci.2009.11.018. [11] Barlat, F., Brem, J.C., Yoon, J.W., Chung, K., Dick, R.E., Lege, D.J., Pourboghrat, F., Choi, S.-H., Chu, E. (2003). Plane stress yield function for aluminum alloy sheets—part 1: theory, Int. J. Plast., 19(9), pp. 1297–1319, DOI: 10.1016/S0749-6419(02)00019-0. [12] Jia, L.-J. (2014). Integration algorithm for a modified Yoshida–Uemori model to simulate cyclic plasticity in extremely large plastic strain ranges up to fracture, Comput. Struct., 145, pp. 36–46, DOI: 10.1016/j.compstruc.2014.08.010. [13] Eggertsen, P.-A., Mattiasson, K. (2011). On the identification of kinematic hardening material parameters for accurate springback predictions, Int. J. Mater. Form., 4(2), pp. 103–120, DOI: 10.1007/s12289-010-1014-7. T R. Fincato et alii, Frattura ed Integrità Strutturale, 57 (2021) 114-126; DOI: 10.3221/IGF-ESIS.57.10 126 [14] Jia, L.-J., Kuwamura, H. (2014). Prediction of Cyclic Behaviors of Mild Steel at Large Plastic Strain Using Coupon Test Results, J. Struct. Eng., 140(2), pp. 04013056, DOI: 10.1061/(ASCE)ST.1943-541X.0000848. [15] Perić, D. (1992). On consistent stress rates in solid mechanics: Computational implications, Int. J. Numer. Methods Eng., 33(4), pp. 799–817, DOI: 10.1002/nme.1620330409. [16] Jiao, Y., Fish, J. (2017). Is an additive decomposition of a rate of deformation and objective stress rates passé?, Comput. Methods Appl. Mech. Eng., 327, pp. 196–225, DOI: 10.1016/j.cma.2017.07.021. [17] Jiao, Y., Fish, J. (2018). On the equivalence between the multiplicative hyper-elasto-plasticity and the additive hypo- elasto-plasticity based on the modified kinetic logarithmic stress rate, Comput. Methods Appl. Mech. Eng., 340, pp. 824–63, DOI: 10.1016/j.cma.2018.06.017. [18] Zbib, H.M., Aifantis, E.C. (1988). On the concept of relative and plastic spins and its implications to large deformation theories. Part I: Hypoelasticity and vertex-type plasticity, Acta Mech., 75(1–4), pp. 15–33, DOI: 10.1007/BF01174625. [19] Tsutsumi, S., Kitamura, T., Fincato, R. (2020). Ductile behaviour of carbon steel for welded structures: Experiments and numerical simulations, J. Constr. Steel Res., 172, DOI: 10.1016/j.jcsr.2020.106185. [20] Gurson, A.L. (1977). Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I—Yield Criteria and Flow Rules for Porous Ductile Media, J. Eng. Mater. Technol., 99(1), pp. 2, DOI: 10.1115/1.3443401. [21] Lemaitre, J. (1985). Coupled elasto-plasticity and damage constitutive equations, Comput. Methods Appl. Mech. Eng., 51(1–3), pp. 31–49, DOI: 10.1016/0045-7825(85)90026-X. [22] de Souza Neto, E.A., Peric, D., Owen, D.R.J. (2008). Computational Methods for Plasticity, 55. [23] Simo, J.C., Hughes, T.J.R. (1998). Computational Inelasticity, vol. 7, New York, Springer-Verlag. [24] Fincato, R., Tsutsumi, S. (2017). Closest-point projection method for the extended subloading surface model, Acta Mech., 228(12), pp. 4213–4233, DOI: 10.1007/s00707-017-1926-0. [25] Fincato, R., Tsutsumi, S. (2018). A numerical study of the return mapping application for the subloading surface model, Eng. Comput., 35(3), pp. 1314–1343, DOI: 10.1108/EC-12-2016-0446. [26] Xie, X., Jiang, W., Chen, J., Zhang, X., Tu, S.-T. (2019). Cyclic hardening/softening behavior of 316L stainless steel at elevated temperature including strain-rate and strain-range dependence: Experimental and damage-coupled constitutive modeling, Int. J. Plast., 114, pp. 196–214, DOI: 10.1016/j.ijplas.2018.11.001.