On finite element implementation of cyclic elastoplasticity: theory, coding, and exemplary problems

In this work the finite element (FE) implementation of the small strain cyclic plasticity is discussed. The family of elastoplastic constitutive models is considered which uses the mixed, kinematic-isotropic hardening rule. It is assumed that the kinematic hardening is governed by the Armstrong–Frederick law. The radial return mapping algorithm is utilized to discretize the general form of the constitutive equation. A relation for the consistent elastoplastic tangent operator is derived. To the best of the author’s knowledge, this formula has not been presented in the literature yet. The obtained set of equations can be used to implement the cyclic plasticity models into numerous commercial or non-commercial FE packages. A user subroutine UMAT (User’s MATerial) has been developed in order to implement the cyclic plasticity model by Yoshida into the open-source FE program CalculiX. The coding is included in the Appendix. It can be easily modified to implement any isotropic hardening rule for which the yield stress is a function of the effective plastic strain. The number of the utilized backstress variables can be easily increased as well. Several validation tests which have been performed in order to verify the code’s performance are discussed.


Introduction
The classical flow plasticity based on the Huber-von Mises-Hencky (HMH) yield criterion is the theory which is the most commonly used for the description of the mechanical properties of metals, e.g., [27]. The constitutive equations of small strain plasticity utilizing the isotropic, the kinematic or the mixed hardening rules are available in practically every finite element (FE) program used for solving engineering problems, e.g., [6,7,9,19]. The standard formulation of the kinematic hardening makes use of the linear rule developed by Prager [20]. This version of the kinematic hardening model is widely used for solving different boundary value problems of elastoplasticity. However, the Prager's model is the most basic one, and often a need arises for a constitutive relation that would describe the material response more accurately. Armstrong and Frederick [2] proposed an elastoplastic constitutive model which utilizes a more elaborate relation to describe the backstress evolution during the kinematic hardening phenomenon. This model is more suitable for the description of material behaviors associated with the cyclic loadings such as the Bauschinger effect. The model by Armstrong and Frederick was further extended by Chaboche and Rousselier by adding the isotropic hardening behavior [4,5]. Thus, a mixed hardening rule was obtained this way. For the purpose of modeling the kinematic hardening, this constitutive model assumes that the total backstress is a sum of several components, each one evolving according to a separate Armstrong-Frederick (A-F) equation. The Voce law [25] was used by Chaboche and Rousselier in order to simulate the isotropic hardening behavior. It was assumed that the total yield stress is a sum of terms which evolve according to separate differential equations of the type defined by the Voce rule. The Chaboche and Rousselier (Ch-R) model and its modifications are often used for simulating the material response in the case of cyclic loadings, e.g., [6,9,19,26]. A special case of the elastoplastic Ch-R model was considered by Yoshida et al. [29] where a single Voce term was used to describe the isotropic hardening behavior, whereas the total backstress associated with the kinematic hardening was assumed as a sum of two components. The evolution of the first backstress component was governed by the A-F law, while the second component was defined by the Prager rule. The rate-independent plasticity Ch-R model was generalized to take into account the viscous effects, cf [4,5]. The Ch-R model is also a basis for the formulation of more elaborate constitutive equations which can take into account such effects as the cyclic shear softening, for instance, e.g., Kowalewski et al. [13].
The Ch-R constitutive model and its various modifications have been implemented into several commercial and non-commercial FE analysis programs, e.g., [6,9,19]. However, usually the details of the FE implementation, such as: the integration of the constitutive equation and the derivation of the tangent operator, are not covered in the available program's documentation. On the other hand, Auricchio et al. [3] discussed the FE implementation of a small strain cyclic plasticity model being a slight modification of the constitutive relation by Chaboche and Rousselier. In the aforementioned work a modified formulation of the HMH yield criterion was utilized. The isotropic hardening was described by a linear function of the effective plastic strain, and a single backstress which evolves according to the A-F rule was used to simulate the kinematic hardening. This approach was further extended by Artioli et al. [1]. Moreover, Kobayashi and Ohno discussed the FE implementation of selected cyclic plasticity models [12]. Kullig and Wippler [14] utilized the implicit backward Euler method and the radial return mapping algorithm to develop an integration scheme for the viscoplastic Chaboche model. A method for calculating the tangent operator was presented as well. The obtained integration algorithm was used to implement the considered constitutive model into the FE system ABAQUS by taking advantage of the UMAT subroutine (User's MATerial), cf [9]. Yang and Feng [28] used the radial return algorithm to implement a viscoplastic material model based on the A-F kinematic hardening law into ABAQUS. A proper UMAT code was written for that purpose.
In this work the FE implementation of the small strain cyclic plasticity is discussed. The radial return mapping algorithm is utilized to develop an integration scheme which is valid for a family of elastoplastic constitutive models assuming the mixed hardening behavior. The general form of the isotropic hardening rule is considered where the yield stress is an arbitrary function of the effective plastic strain. It is assumed that the kinematic hardening behavior is determined by a number of backstress variables which evolve according to the separate A-F equations. In the special case, when the function describing the isotropic hardening behavior is assumed in the form of the Voce law, the considered constitutive model takes the form of the Ch-R model.
A general relation for the consistent fourth-order tangent operator is derived. To the best of the author's knowledge, this formula for the tangent operator of the considered model formulation has not been presented in the literature before. The fact that the fourth-order operator is consistent, i.e., in agreement with the radial return mapping algorithm used for the integration of the model equations, guarantees a quadratic rate of convergence during the numerical computations when the operator is entirely coded in an FE program, cf e.g., [1].
The developed model integration algorithm along with the tangent operator is utilized to implement the elastoplastic model proposed by Yoshida et al. [29] into the FE program CalculiX. For that purpose, the user subroutine UMAT has been written. Numerous FE simulations have been performed in order to verify the performance of the developed code. The obtained results are presented in this work. The UMAT code has been attached in the Appendix Section. It has a general form and can be easily modified to incorporate some different isotropic hardening laws (Ludwik [16], Swift [24] or numerous Voce terms etc.). The number of programmed backstresses defined by the A-F or the Prager laws can be easily increased as well.

Basic notions
The total stress tensor σ can be decomposed into the volumetric stress p (opposite to the pressure) and the deviatoric stress s, e.g., [23], i.e., The volumetric stress and the stress deviator are associated with the elastic strain by the following relationships: where K is the bulk (Helmholtz) modulus and μ = G is the shear (Kirchhoff) modulus, whereas e and e e are the elastic dilatation and the elastic deviatoric strain, respectively, i.e., e = tr ε e , e e = ε e − 1 3 where ε e is the total elastic strain tensor. The shear and the bulk moduli are related to the engineering constants, i.e., Young's modulus E and Poisson's ratio ν, via the following equalities: It is assumed that the total small strain tensor ε can be additively decomposed into its elastic and plastic components, i.e., ε e and ε p , respectively [20,23,27]: where ε p is the plastic strain tensor. The volume change is treated as fully elastic, thus tr ε p = 0, so that ε p = e p . It follows that in Eqs. (2) e = = tr ε while e e = e − e p . Thus, another equation which governs the evolution of the plastic strain tensor is required. The Huber-von Mises-Hencky (HMH) yield criterion has to be satisfied when plastic flow occurs [4,5]: where X is the backstress tensor, k is the initial yield stress, while R(ē p ) is a function responsible for the isotropic hardening behavior (the expansion or the reduction in size of the yield function F(σ , X,ē p )). The effective plastic strainē p is given by the following formulas: The evolution of the plastic strain is governed by the equation: where dλ is the plastic multiplier. After some manipulation, it can be found from Eq. (8) that dλ = dē p . Thus: wheren is the normalized effective stress tensor. The total backstress tensor X is assumed to be a sum of M components X (i) (i = 1, 2, . . . , M) [4,5], i.e., which evolve according to the separate Armstrong-Frederick (A-F) equations, cf [2]: where C i and γ i are the kinematic hardening parameters. Some specific forms of R(ē p ) are considered further in the text. In the following Section the discretization of the constitutive model equations is discussed.

Numerical integration of the constitutive model
The radial return mapping algorithm is utilized for the integration of the constitutive model, e.g., [1,3,14,22]. The deviatoric stress in the increment n + 1 is given by the following relation: thus: s n+1 = s pr n+1 − 2μ e p , s pr n+1 = s n + 2μ e, (13.1,2) with s pr n+1 being the predictor stress. Equations (11) are discretized by replacing the differentials with the finite differences, i.e., After some rearrangements Eq. (14) can be rewritten as follows: Equations (9) can be written in incremental form, i.e., After some manipulation Eq. (17) can be written as: Using Eq. (16.2) in Eq. (18) and some rearrangements lead to: where when the material is yielding. The following notations are introduced: and After inserting Eqs. (20)- (22) into Eq. (19) it is found that: By performing some manipulations Eq. (23) can be converted into a scalar equation, i.e., where J 2 (Z) = 3 2 Z · Z andē p n+1 =ē p n + ē p . Equation (24) is a nonlinear algebraic equation which has to be solved numerically for the effective plastic strain increment ē p . The solution of Eq. (24) by means of Newton's iterative method is discussed in the Appendix. After determining the value of ē p the backstress components can be updated according to Eq. (15.1) and the total backstress is calculated using Eq. (10). Subsequently, the deviatoric stress is updated using Eq. (23), i.e., It is important to note that using Eqs. (20), (23), and (24) it can be found that:

Consistent tangent operator
Taking the variation of Eq. (13.1) with respect to all quantities gives: where the formulas for the variations δē p and δn n+1 have to be found. For the purpose of finding δē p , the variation of Eq. (24) is taken, i.e., with δ R(ē p ) = d R(ē p ) dē p δē p which depends on the specific form of R(ē p ) (see Table 1), whereas δα( ē p ) = dα( ē p ) d ē p δē p takes the form: see Appendix for the derivation. The variation δ J 2 (Z) is given by the following equation (see Appendix): Substituting Eqs. (29) and (30) into Eq. (28) and solving it for δē p leads to the following relation: It is convenient to introduce the so-called effective hardening modulus: Thus, using Eq. (32) the formula for the effective plastic strain variation can be written as: The variation of the normalized effective stress is given as (see Appendix): which after some rearrangements can be expressed as: The effective shear modulus μ * is introduced as: Moreover, the following relations are utilized: and δs Substituting Eqs. (37)-(39) into Eq. (36) yields: where the fact that 1 ·n n+1 = tr(n n+1 ) = 0 has been utilized. The variation of the total stress tensor at the increment n + 1, i.e., σ n+1 , is given as: A relation between the stress and the strain variations can be obtained by inserting Eqs. (40) and (41.2) into Eq. (41.1), i.e., where C C C e− p n+1 is the consistent algorithmic tangent operator which is given by the following relation: being the effective Lamé parameter. It should be emphasized that the tangent operator given by Eq. (43) is non-symmetric, cf [3,14]. For M = 1 and γ 1 = 0 the tangent operator simplifies to the form corresponding to the mixed hardening model (with the kinematic hardening governed by the rule by Prager, cf [11]). Moreover, if it is further assumed that R(ē p ) = 0 the tangent operator of the Prager's kinematic hardening model is obtained, cf [22]. If one assumes M = 0 in Eq. (43) the tangent operator is reduced to the one of the isotropic hardening model [22]. If the equivalent HMH predictor stress does not satisfy the condition , the tangent operator takes the form of the elasticity tensor, i.e., where λ = K − 2 3 μ is the Lamé constant.

User material (UMAT) subroutine for cyclic elastoplasticity
Most of the commercial and non-commercial FE programs offer the option of using a user-written subroutine to define a constitutive relation which is unavailable in the program's material library, e.g., [6,7,9,19]. Usually such subroutines require from the user to define the stress update rule and the material tangent operator also known as the material Jacobian. The FE program CalculiX utilizes alternatively two different interfaces which can be used for coding a user-defined constitutive law, cf [7]. The first is the ABAQUS interface with the second being the CalculiX native interface. In this study the stress update algorithm described in Section 3 and the tangent operator given by Eq. (43) were used to implement the small strain cyclic plasticity constitutive model into the FE program CalculiX via the UMAT subroutine (User MATerial). It was decided to utilize the CalculiX native interface due to its better performance compared to the ABAQUS UMAT used under CalculiX. The subroutine UMAT uses the small mechanical strain tensor components (the "emec" column matrix) as an input which is further used to calculate the stress tensor component matrix ("stre") and the components of the tangent operator stored in a column matrix ("stiff"). It is assumed that the tangent operator is symmetric [7], i.e., when the tangent operator is written as a 6 × 6 matrix only the upper half of the components should be defined. The non-symmetric tangent operator tensors should be symmetrized as shown in the exemplary source files provided with CalculiX. It has been demonstrated in the literature that the performed symmetrization has no negative influence on the accuracy of the FE computations and very limited influence on their robustness, e.g., [10,12]. Below the indexes of the fourth-order tensor components have been listed in the 6 × 6 matrix (the indexes of the "stiff" column matrix components are given in the brackets; the Voigt notation is not used): The subroutine UMAT calculates the components of stress and material tangent operator at each Gauss integration point. These quantities are subsequently used by CalculiX to form up the element stiffness matrix. Finally, the global stiffness matrix is assembled by CalculiX using the element stiffness matrices. This procedure is repeated during every iteration of the Newton-Raphson process for all increments of the analysis (Fig. 1).
The UMAT subroutine is written in Fortran 77 language. In order to use it with CalculiX the code should be compiled to the dynamic link library file (DLL) 1  where "YOSHIDA" is the name of the DLL file containing the compiled UMAT code for the elastoplastic model by Yoshida [29]. This UMAT code can be found in the Appendix Section. The calculations performed by the developed subroutine follow the list in the box below. In the case of the Yoshida elastoplastic model M = 2 with γ 2 = 0, whereas the isotropic hardening behavior is simulated using the rule by Voce, cf [29].
Marciniak and Kuczyński [17] Aē p +B The number of the utilized backstresses can be easily increased by a proper modification to the attached code. The UMAT subroutine has been written in a general form so that any other function than Voce defining the isotropic hardening behavior can be easily implemented by modifying the code. Some alternative forms of the function R(ē p ) and its derivative are listed in Table 1.
It should be emphasized that the maximum number of material parameters which can be listed in a single row in the *Material card is eight. If there are more parameters, they should be listed in the subsequent rows. Moreover, if the number of the model's parameters is precisely eight (as in the example given above) the temperature should be given in the subsequent row of the input file, regardless of the type of analysis [7].

Exemplary problems
Below some exemplary simulations are discussed which were utilized to verify the performance of the developed UMAT code. The verification tests included simulating processes with homogenous strain and stress fields: uniaxial tension/compression (UT/UC), equibiaxial tension/compression (BT/BC) and simple shear (SS). The radial return mapping algorithm was used to derive equation sets which describe the aforementioned processes. These equations were utilized to write simulation programs in Scilab software [21] which were further used to verify the results of the FE simulations performed with CalculiX. In addition to the problems mentioned above some more complex examples were considered as well.

Uniaxial tension/compression (UT/UC)
A uniaxial stress state is assumed. The matrices containing the components of the total stress tensor, its deviator (according to Eq. (1.3)) and the backstress tensor are given as: Thus: and the HMH equivalent stress takes the form: It follows from the generalized Hooke's law given by Eqs. (1.1) and (2) that the axial elastic strain component is where according to the assumed boundary conditions: thus: The axial stress component in the increment n + 1 is given as: with the axial elastic strain component The total strain and plastic strain matrices take the form: where ε a is the total axial strain, ε l is the total lateral strain, whereas ε p is the axial plastic strain component. It follows from Eqs. (54) and (55) that Eq. (53) can be rewritten as: The components of the Z stress are assumed in the following form: Thus, it follows from Eq. (26) that the axial component of the normalized effective stress is: After substituting Eq. (58) into Eq. (16.1) it is found that: Inserting Eq. (59) into Eq. (56.1) results in According to Eq. (15.1) the axial component of the i-th backstress in the increment n + 1 is while it follows from Eq. (10) that the total axial backstress is whereas X 11 = 2 3 X and X (i) After inserting Eq. (59) into Eq. (61) it is found that Subtracting Eq. (63) from Eq. (60) and substituting Eq. (64) gives: It follows from Eqs. (26), (47), and (57) that According to Eq. (20) the following equality should hold when the material is actively yielding: Substituting Eq. (68) into Eq. (67) leads to: Thus, after inserting Eq. (69) into Eq. (65) and some rearrangements it is found that: whereē p n+1 =ē p n + ē p . The following notation is used: It follows from Eqs. (22), (47), and (57) that Substituting Eqs. (71) and (72) into Eq. (70) after some rearrangements gives: The following equation is also valid: Inserting Eq. (68) into Eq. (74) and further rearrangement yields: which is the nonlinear algebraic equation that should be solved for the effective plastic strain increment ē p . If the isotropic hardening function is assumed in the form proposed by Chaboche and Rousselier, then: The derived set of equations was utilized to write a Scilab program designated for simulating uniaxial tension and compression processes. The computations performed by the program have been listed in the box below. The Scilab function fsolve [21] was used for solving the nonlinear algebraic equation given by Eq.
(75). The developed program was used to verify the results of the FE simulations performed using the UMAT code.
The data for DP1000 steel by Zimniak and Wiewiórska [30] were utilized to determine the material parameters of the elastoplastic constitutive model by Yoshida. The evaluated parameter values are listed in Table 2. This material parameter set was used for the numerical simulations both in CalculiX and Scilab. In the case of the Yoshida model M = 2 in Eq. (62) with γ 2 = 0 and N = 1 in Eq. (76).
The simulation of the uniaxial tension/compression process was performed in CalculiX using a cubic geometrical model with the dimensions 1mm×1mm×1mm. The cube was meshed with a single C3D8 element 2 . In Fig. 2a the applied boundary conditions have been illustrated. A kinematic excitation was assumed. A ramp displacement δ in the direction "1" of the rectangular coordinate system was applied to the cube's frontal face ABC D (see Fig. 2c). The following boundary conditions were applied to the other faces of the cube: a zero displacement in the direction "1" was set on the face E FG H, a zero displacement in the direction "2" on the face AE H D, and a zero displacement in the direction "3" on the face AB F E. After reaching the maximum value of δ = 0.07 mm the displacement started decreasing linearly with the analysis time. The simulation ended with δ = 0 mm. The axial strain is calculated as: ε a = δ/1 = δ.
(d) Update strains:    The results obtained using the finite element method (FEM) in CalculiX were compared to those generated by the simulation using the Scilab program. In Fig. 3a the comparison of the stress-axial strain responses produced independently by the two aforementioned methods can be seen. In Fig. 3b the obtained plots of the axial strain (ε a ) vs the lateral strain (ε l ) are compared. An excellent agreement is found in both cases.
Another uniaxial tension/compression simulation was performed using the same FE model and set of boundary conditions with a more complex loading history. In this approach the displacement of the cube's ABC D face was defined using the increasing saw function. It can be seen in Fig. 4a and b that again an excellent agreement was found between the results produced by CalculiX and Scilab. In Fig. 4c the assumed complex loading history can be seen.
Both of the described FE simulations were repeated for the cube meshed with twenty-seven C3D8 elements. Again a very good agreement was observed between the results generated using CalculiX and Scilab. The simulations described above were also performed for different types of finite elements, i.e., C3D20 3 , C3D4 4 and C3D10 5 . No decrease in the performance of the developed UMAT code was observed.

Equibiaxial tension/compression (BT/BC)
An equibiaxial stress state is assumed. The stress tensor components and the volumetric stress are given as: 3 Cubic, three-dimensional, twenty nodes. 4 Tetrahedral, three-dimensional, four nodes. 5 Tetrahedral, three-dimensional, ten nodes.
It follows that the components of the stress deviator take the form: The backstress and the auxiliary stress Z are assumed in a similar form to the stress deviator, i.e., whereas, due to the plastic incompressibility (tr(ε p ) = 0), the plastic strain has the following components: It follows from Eqs. (78) and (79.1) that thus the HMH equivalent stress takes the form: According to the generalized Hooke's law as given by Eqs. (1.1) and (2) after some rearrangements the elastic strain components are given as: After taking into account that σ 11 = σ 22 = σ and σ 33 = 0 and some further manipulations it is found that: The matrix of the components of the total strain tensor takes the form: It follows from Eqs. (85) and (5) that: Thus, Eq. (84.1) can be written in the following incremental form: so that the stress in the increment n + 1 is given as: After introducing the predictor stress Eq. (88) takes the form: It follows from Eqs. (79.1,2) and (26) that According to Eqs. (16.1) and (90.2) the plastic strain increment ε p 11 is given by the following relation: Assuming X It follows from Eqs.
The Z 11 component of the auxiliary stress is given as: Thus, according to Eqs. (79.2) and (96): After substituting Eqs. (94), (95) and (97.2) into Eq. (93) and some manipulations it is found that: where the following notation is adapted: with w i (i = 1, 2, . . . , M) given by Eq. (15.2). Equation (98) can be transformed into the form: which is the nonlinear algebraic equation that has to be solved numerically for the effective plastic strain increment ē p . The following relation which follows from Eq. (98) is used for updating the stress when the plastic strain increment is determined: In the case of the Yoshida model M = 2 with γ 2 = 0, whereas the isotropic hardening is described using the model by Voce (N = 1 in Eq. (76)). The derived equations were used to write a Scilab program designated for simulating the equibiaxial tension/compression processes. The calculations performed by the program during every computational step have been gathered in the box below. The fsolve function offered by Scilab was utilized for solving Eq. (100). The developed Scilab program was used to verify the results of the FE simulations performed using the UMAT code. In Table 2 the material parameter values which were used for the simulations are gathered.
The simulation of the equibiaxial tension/compression process was performed in CalculiX using a cubic geometrical model with the dimensions 1mm×1mm×1mm. The cube was meshed with a single C3D8 element. In Fig. 2b the applied boundary conditions have been illustrated. A kinematic excitation was assumed. A ramp displacement δ in the direction "1" of the rectangular coordinate system was applied to the two of the cube's faces, i.e., ABC D and B FGC (see Fig. 2c). The following boundary conditions were applied to the other faces of the cube: a zero displacement in the direction "1" was set on the face E FG H, a zero displacement in the direction "2" on the face AE H D, and a zero displacement in the direction "3" on the face AB F E. After reaching the maximum value of δ = 0.05 mm the displacement started decreasing linearly with the analysis time. The simulation ended when δ = 0 mm.
The results obtained in CalculiX were compared with those generated by Scilab program. In Fig. 5a the stress-axial strain (the strain in the axes of tension) response produced by FEM simulation was compared with that generated by the Scilab program. In Fig. 5b the obtained plots of the lateral strain vs the axial strain are compared. In both cases an excellent agreement is found between the results generated by CalculiX and Scilab. The same analysis was repeated for a larger number and different types of finite elements (C3D4, C3D10, C3D20). Again, a very good agreement between the results obtained using CalculiX and Scilab was found.

Simple shear (SS)
In the case of SS process the stress, the backstress, and the auxiliary stress tensors have the following components: It follows from Eqs. (2) and (5) that the total shear stress in the n + 1 increment can be expressed as where It follows from Eq. (113) that Inserting Eq. (112) into Eq. (115) gives which is the nonlinear algebraic equation that has to be solved numerically for ē p at every computational step of the simulation. When the law by Chaboche and Rousselier is applied to simulate the isotropic hardening behavior, the function R(ē p n + ē p ) is given by Eq. (76). Substituting Eq. (108) into Eq. (104.2) results in the following formula for updating the backstress: The derived set of equations was utilized to write a Scilab program which can be used for performing simulations of the SS processes. The calculations performed by the Scilab program at every computational step of the numerical simulation are listed in the box. The fsolve function was utilized for solving Eq. (116). Again, in the considered case of the Yoshida model M = 2, γ 2 = 0, while the isotropic hardening behavior is governed by the Voce rule (N = 1 in Eq. (76)).
The developed Scilab program was used to validate the results generated by CalculiX during the SS process simulation. The material parameter set given in Table 2 was utilized. The boundary conditions used for the simulation performed in CalculiX are illustrated in Fig. 6a. A cube was meshed with a single C3D8 finite element. The displacements on the ABC D face of the cube (see Fig. 6b) were set to zero. The displacements in the directions "2" and "3" were set to zero on the face E FG H. The displacement δ of the face E FG H in the direction "1" was used as the kinematic excitation which increased linearly up to the maximum value of δ = 0.2 mm. After reaching the maximum the displacement value started decreasing linearly to zero. Since δ = tan γ ≈ γ (Fig. 6), it follows that the shear strain ε = γ /2 = δ/2. In Fig. 7 the comparison of the results produced by Scilab and CalculiX is presented. An excellent agreement was found. The simulation was repeated for the cube meshed using C3D4, C3D10, and C3D20 elements. In all the considered cases a very good agreement was found between the results generated by the Scilab program and those obtained from the FE simulation in CalculiX.  In order to investigate the accuracy of the FE computations performed using the developed UMAT subroutine a simulation of cyclic biaxial tension/compression was performed in CalculiX. Again, a cube with the edge length of 1mm, meshed with a single C3D8 element, was used for the simulation. A kinematic excitation was assumed. A displacement δ 1 in the direction "1" of the rectangular coordinate system was applied to the cube's face ABC D, Fig. 2c. The time history of δ 1 is shown in Fig. 8. Another displacement δ 2 in the direction "2" was defined on the cube's B FGC face. The history of δ 2 can be seen in Fig. 8. The following boundary conditions were applied to the other faces of the cube: a zero displacement in the direction "1" was set on the face E FG H, a zero displacement in the direction "2" on the face AE H D, and a zero displacement in the direction "3" on the face AB F E.
The simulation was repeated for different sizes of the fixed time increment t, i.e., 1E-2, 1E-3, 1E-4, and 1E-5 s. The FE solution obtained for t = 1E − 5 s was assumed as exact. The following normalized error function was used to estimate the accuracy of FE computations for different incrementation size: with σ exa being the Cauchy stress tensor calculated for t = 1E − 5 s, whereas σ com is the stress tensor computed for another increment size for which the error value is being calculated. In Fig. 9a and b the stressstrain histories that were computed for the selected increment sizes are shown. In Fig. 9c Fig. 10a. A ramp displacement δ of the bar's upper face was used as a kinematic excitation. The bar was elongated until the maximum value of δ = 1 mm was reached. The bar was meshed using C3D8 elements (see Fig. 10b). The material parameter values gathered in Table 2 were used for the Yoshida elastoplastic model defined by subroutine UMAT. In Fig. 11 the results of the FE simulation are presented. The largest residual force (LRF) values recorded during the subsequent iterations are gathered in Table 3. It can be seen that for some increments the quadratic convergence of solution was achieved.

Cyclic loading
In the second approach the same set of boundary conditions as before was used; however, the time history of the displacement δ (Fig. 10a) was defined using a triangular wave function, Fig. 12a. The simulation was performed for different fixed increment sizes, i.e., 5E-3, 2.5E-3, and 1E-3 s. The total reaction force F at the bar's upper face was calculated for each increment. The plot of F vs. δ for selected incrementations can be seen in Fig. 12b. In Table 4 the analysis time values obtained for different increment sizes are gathered 6 . It is seen that the analysis time is approximately inversely proportional to the increment size.

Hollow cylinder under cyclic tension/compression and torsion
A hollow cylinder subjected to cyclic axial tension/compression and torsion was considered. The cylinder was meshed with 19440 C3D8 elements, Fig. 13a. One end of the cylinder was assumed fixed with all displacements set to zero. An axial displacement was defined on the cylinder's frontal face. Moreover, a rotation of the frontal face around the cylinder's axis was applied. The assumed time histories of the axial displacement δ and the angle of twist α are shown in Fig. 13b. The maximum values were δ max = 1.08 mm for the axial displacement and α max = 9.8 o for the angle of twist.
The simulation was performed in two variants using the automatic incrementation option. For the first variant the maximum time increment size t max was set to 1E-1 s, while for the second variant it was set to 2.5E-2 s. The computed results that were saved for the selected finite element no. 29983 (Fig. 13a) can be seen in Fig. 14. Very similar results were obtained for both considered incrementation methods.

Conclusions
In this work the FE implementation of cyclic elastoplasticity was discussed. The radial return mapping algorithm was used to develop a numerical integration scheme for a group of constitutive equations based on the     [4,5]. A fourth-order tangent operator was derived which is consistent with the utilized algorithm of integrating the constitutive model. To the best of the author's knowledge this particular form of the tangent operator has not been presented in the literature before.
A subroutine UMAT was developed which allows to implement the elastoplastic model by Yoshida [29] into the FE program CalculiX. The UMAT code has a general form which allows to easily modify the isotropic hardening rule used by the model. The number of the utilized backstress variables can be increased easily as well. Many validation tests were performed in order to verify the performance of the developed UMAT subroutine. An excellent agreement was found between the results produced by UMAT and those of the verification programs. What is more, a quadratic convergence was observed during the Newton-Raphson iterative process.