A new method for identification of cyclic plasticity model parameters

In this study, a new method for determining the material parameters of cyclic plasticity is presented. The method can be applied to evaluate the model parameters from any loading histories measured experimentally. The experimental data require basic processing only to be utilized. The method can be applied to calibrate the parameters of different elastoplastic models such as the Chaboche–Rousselier (Ch–R) constitutive equation or other model formulations which use different rules of isotropic hardening. The developed method was utilized to evaluate the material parameters of copper for a selected group of constitutive models. It is shown that among the considered model formulations a very good description of the mechanical properties of copper is achieved for the Ch–R model with two Voce terms used for simulating the isotropic hardening and two backstress variables utilized for capturing the kinematic hardening behavior. Furthermore, it is demonstrated that a model calibrated using the cyclic tension/compression data is able to properly capture the material response in torsion. Similarly, when the constitutive parameters are determined using the cyclic torsion data the model is able to properly reproduce the material behavior in tension/compression. It is concluded that for the considered type of constitutive equations the material parameters can be identified from a single mechanical test. The proposed methodology was validated using the relations derived analytically


Introduction
Despite the continuous development of new functional materials, such as polymers and composites, e.g., [1,2], metals are still widely used in different branches of industry. A case where metallic structural elements are subjected to cyclic loadings is common. Thus, a proper description of material response which takes into account such behaviors as the Bauschinger effect, for instance, is indispensable.
Over the years, many constitutive equations were developed with the purpose of describing mechanical properties of metals in both the elastic and the elastoplastic ranges of strains. The constitutive models which are based on the notions of yield condition and flow rule are the most commonly used. Usually, the Huber-von Mises-Hencky (HMH) yield criterion or its modifications are utilized to simulate the material's plastic flow and the hardening behavior [3]. A constitutive equation which assumes isotropic hardening of the material under loading is one of the simplest model formulations. The yield stress is defined as a function of the effective plastic strain. Usually, a piece-wise linear relation between the yield stress and the effective plastic strain is defined. However, sometimes nonlinear relations are used such as the one proposed by Voce [4]. Another classical model of the flow plasticity theory utilizes the so-called kinematic hardening behavior. In the simplest case the hardening rule developed by Prager [5] is employed for that purpose. This model formulation allows to take into account the Bauschinger effect.
In most cases, the simple classical models of the flow theory of plasticity fail to properly capture material response when it is subjected to more complex and cyclic loadings. Thus, there is a need for more sophisticated constitutive models which are able to correctly describe the material behavior in these conditions. Armstrong and Frederick [6] developed a more advanced kinematic hardening rule which describes the evolution of the so-called backstress variable. This model formulation was subsequently generalized by Chaboche and Rousselier [7,8] whose model utilizes the Armstrong-Frederick (A-F) equation to capture the kinematic hardening phenomenon and the Voce's rule to take into account the isotropic hardening (or softening) effect. Multiple backstress variables and Voce exponential terms can be utilized to properly describe the aforementioned phenomena. A modification of the Chaboche-Rousselier (Ch-R) model was used by Yoshida et al. [9]. Moreover, a generalization of the Ch-R model was proposed by Chaboche to include the viscous effects [10].
The determination of constitutive model's material parameters presents itself a very important problem. The parameters need to be properly evaluated for a specific material under study, so that, the model will be able to produce correct results of the analytical and numerical engineering calculations.
Several analytical or quasi-analytical methods exist which allow to estimate the values of the Ch-R model parameters based on the experimental results of the loading-unloading processes. On the other hand, a number of parameter identification algorithms were proposed which are based on the least-squares method. Some of the proposed parameter determination techniques utilize both the analytically derived relationships and the least squares method, e.g., [11][12][13]. Zimniak and Wiewiórska [14] developed a parameter determination method for the Ch-R model in which a plot of the backstress variable as a function of the effective plastic strain is generated from the cyclic tension/compression data. Subsequently, the experimental measurements of the backstress are fitted using an analytically derived relation. In the following step, the remaining parameters responsible for the isotropic hardening are determined by approximating the curve of the yield stress vs. plastic strain. Wójcik and Skrzat [15,16] developed a material parameter evaluation algorithm which utilizes some analytical formulas for calculating the theoretical stress values that are further used within the least squares optimization process. In this concept the fuzzy logic methods are applied in order to improve the curve-fitting quality.
A common approach in the determination of Ch-R model parameters is performing a finite element (FE) simulation in order to calculate the theoretical stress values which are further used in the least-squares method. This approach is usually coupled with usage of a genetic algorithm for the optimization of material parameter values [17][18][19].
In this work, a new method is proposed for the determination of elastoplastic parameters of the Chaboche-Rousselier constitutive model. The searched parameters are evaluated by minimizing the total square error. The main novel idea of the presented approach is calculating the theoretical stress values numerically using the radial return mapping algorithm. The presented method has the following advantages: • It reduces the amount of experimental data processing which is necessary for determination of material parameter values. • The parameters are determined fast in a single identification stage. • The method can be used for evaluating the material parameters from any experimentally registered loading histories. • It can be applied to the modifications of Ch-R model that use different isotropic hardening rules.
Correctness of the conducted numerical calculations has been validated by comparing the obtained results to the predictions of the derived analytical equations describing tension/compression and torsion processes. A number of cyclic tension/compression and torsion/ reverse torsion tests were conducted on copper specimens. The experimental measurements were used to identify the material parameters of copper by taking advantage of the newly developed method. It is concluded that very good results are obtained when the mechanical properties of copper are simulated with a model employing two Voce terms to capture the isotropic hardening and two backstress variables to describe the kinematic hardening behavior.
Furthermore, it is shown that a model calibrated using the cyclic tension/compression data is able to accurately describe the material response in torsion. Similarly, the material parameters identified using the cyclic torsion data allow to correctly predict the material behavior in tension/ compression. Thus, it is concluded that for the considered group of constitutive equations the material parameters can be evaluated based on a single mechanical test (tension/compression or torsion).

Basic notions
The total stress tensor is a sum of the volumetric stress p and the stress deviator , i.e., where the elastic volume change e = tr ( e ) with tr (•) being the trace operator, e is the elastic strain deviator, whereas K and are the bulk and shear moduli, respectively.
In the classical theory of plasticity the total strain tensor is taken to be the sum of the elastic strain and the plastic strain, e.g., [3], i.e., An additional relationship which determines the evolution of the plastic strain is required. The Huber-von Mises-Hencky (HMH) yield condition is utilized in the following form: where is the backstress variable which governs the kinematic hardening behavior, k is the initial yield stress, R is a function specifying the isotropic hardening, whereas ē p is the effective plastic strain that is given as where ̄ is the normalized effective stress. In the case of Chaboche-Rousselier (Ch-R) model the function R is assumed as a sum of terms in the form proposed by Voce [4,7,8]: where Q i is the saturated value of R i , whereas b i governs the rate at which the saturation of R i is reached ( i = 1, 2, … , N).
The total backstress tensor can be viewed as a residual stress on the material's microstructural level. It is assumed as a sum of components (i) ( i = 1, 2, … , M ), i.e., which evolve according to separate equations of the form proposed by Armstrong and Frederick [6]: where C i and i are the material parameters responsible for strain hardening and dynamic recovery, respectively. The ratio C i ∕ i is equal to the saturated, steady-state value of the axial component of (i) in the case of uniaxial tension. Below the algorithm used for the integration of the constitutive model is discussed.

Model discretization
The radial return mapping algorithm is used to integrate the considered elastoplastic constitutive equation, cf [20,21]. At every analysis step the so-called predictor stress is calculated. For the increment no. n + 1 ) is satisfied, then the stress deviator in the increment n + 1 is taken as equal to the predictor stress, i.e., n+1 = pr n+1 . If the aforementioned condition is not satisfied, the plastic flow occurs and the effective plastic strain increment Δē p needs to be determined by solving the following nonlinear algebraic equation: where the function (Δē p ) is given as with and After calculating the effective plastic strain increment the strain and the backstress variables are updated according to the following relations: Subsequently, the total deviatoric stress is updated, i.e., The derivation of specific equation sets for the cases of uniaxial tension/compression and torsion/reverse torsion processes can be found in [21]. These equations were utilized for the purpose of approximating the experimental data and are summarized in the following sections.

Material, specimen and testing machine
Experiments were performed on specimens made of pure copper M1E (notation according to the Polish Standard , respectively. Two separate servo-controller units connected to the computer of the INSTRON loading system can independently apply controlled axial load and torsional moment. The hydraulic pressure in the actuators comes from two servo-valves operated by servocontrollers provided with set-point control signals from the computer. A multiple analogue-to-digital converter (TestStar II) feeds the computer with the signals of: axial displacement machine piston, rotation of the grip fixture of a specimen, axial force, twisting moment, axial strain and torsional strain.
The axial force and the torque applied to the specimen were measured using load cells incorporated in the machine. The software developed for these tests enabled the maintenance of constant strain rates during plastic loading and the resulting stress-strain responses were recorded by an acquisition unit. It was connected to the control computer enabling both direct on-line observations of the experimental results and also their saving onto the hard disk of the computer during each test.
The strains during standard tensile tests were captured by means of the MTS extensometer of the gauge length equal to 25 [mm]. The strains during low cyclic fatigue (LCF) tests were measured with the use of two temperature compensated 45 • rosette strain gauges (type: EP-08-125RA-120 with gauge factor of 2.06 ± 1 [%] , temperature range from [ −75 • C] to +205 [ • C], and grid resistance of 120.0 [ Ω ] ±0.2 [%] ) bonded to the outer surface of the specimen located on the opposite sides of the specimen gauge length and two additional separate strain gauges located on a specially designed semi-ring, which were used for temperature compensation in axial direction. The three stacked strain gauges in the 45 • rosette were arranged in such a manner that one strain gauge was aligned with the longitudinal axis of the specimen while the other two gauges made a 45 • angle, symmetric with respect to the longitudinal gauge. Thus, such a strain measurement system enables independent monitoring of axial and shear strains by means of two full bridge circuits. These strain measurement circuits were connected to the INSTRON measurement systems. At all tests the hoop strains were also registered by means of the additional half bridge circuit of two strain gauges bonded to the specimen and connected to the INSTRON measurement system. Before running the tests, all circuits were calibrated using a highly sensitive Hottinger tensometric bridge (UGR 60). All tests were controlled by a TestStar II using such software as TestWare-SX v. 4.0D.
To fix the thin-walled tubular specimen in the testing machine a special gripping system was designed. It enabled elimination of the backlashes in assembly of the specimen and its automatic alignment during loading.

Program of tests
The comprehensive research program for strength and LCF testing included: 1. Static tensile tests (3 tests) for pure copper in the asreceived state carried out to determine basic mechanical parameters; Fig. 1 Thin-walled tubular specimen 2. LCF tests under symmetric tension-compression and torsion-reverse torsion in order to determine parameters necessary for cyclic models evaluation.
Uniaxial tensile tests were performed at strain rate equal to 2 × 10 −4 s −1 . Details of the experimental setup are shown in Fig. 2. The range of fatigue loads was established on the basis of the yield strength R 0.2 determined from the uniaxial tensile test. LCF tests were carried out according to a strain controlled loading program under frequency of cycles equal to 0.02 Hz and equivalent strain of 0.8% . For both types of cyclic loading (tension-compression and torsion-reverse torsion) the number of cycles was the same (15 full cycles).

Results of tests
A range of mechanical parameters determined from tensile characteristic of copper are presented in Table 1.
The experimental results of LCF tests clearly demonstrated, that the copper exhibited a significant softening effect that is similar for both types of cyclic loading. The hysteresis loops obtained for both programs will be discussed wider during presentation of comparison of the experimental data to predictions achieved by the models proposed.

Approximation of cyclic tension/ compression data
Equations (10-16) that were derived using the radial return mapping algorithm were subsequently transformed into a specific form which is correct for the uniaxial tension/compression process. The obtained equation set is presented below. It was utilized for calculating the theoretical stress values during the least squares approximation of experimental data.

Calculating theoretical stress values
In the case of uniaxial tension/compression process the stress state has one non-zero axial component, i.e., 11 = while 22 = 33 = 0 , which leads to the stress deviator and the backstress as given below: At every increment, the predictor stress is calculated first. Thus, for the computational step n + 1: with Δ a being the increment of the axial strain. Subsequently, the yield criterion given by Eq. (3)

Identification of material parameters
The experimental measurements gathered in the cyclic tension/compression test along with the derived theoretical relationships were utilized for the determination of material parameters of Chaboche-Rousselier (Ch-R) model. For that purpose, the least squares method was used. The minimization of the total square error was performed using the Scilab software and the fminsearch function which is available in this program [22]. The following error function to be minimized was defined: where j is the number of a collocation point, n is the total number of used collocation points, is the theoretical stress value calculated with Eqs. (17)(18)(19)(20)(21)(22), ̃ is the experimentally measured stress and is the matrix of material parameter values being optimized. The previously determined value of the Young's modulus, i.e. E = 113 [GPa] (cf Table 1) was used in the calculations, whereas the initial yield stress k = 145 [MPa] was determined from the uniaxial test data. The Poisson's ratio = 0.32 [-] was estimated based on the values given in the available databases [23]. The determined values of material constants were used during the process of evaluating the remaining parameters of the Ch-R model.
The material parameters were identified for three different versions of the Ch-R constitutive model, i.e.:  Fig. 3 Comparison of experimental results obtained for copper (cyclic tension/compression) and theoretical predictions of Chaboche-Rousselier elastoplastic model with one isotropic hardening term ( N = 1 ) and one backstress variable ( M = 1) • A modification of the model described above obtained by adding an additional backstress variable which follows the Prager's rule (this model formulation was used by Yoshida et al. [9]); • A model with two Voce terms ( N = 2 ) and two backstress variables ( M = 2 ) that evolve according to two separate A-F equations; • A model with three Voce terms ( N = 3 ) and three backstress variables ( M = 3 ) that evolve according to three separate A-F equations.
The curve-fitting results obtained for the first version of the model ( N = 1 , M = 1 ) can be seen in Fig 3. The material parameters of the Ch-R model are gathered in Table 2.
Only slightly better approximation was achieved for the model formulation that was proposed by Yoshida. Thus, the plot showing the approximation is not included here, whereas the identified material parameter values are gathered in Table 3.
A significant improvement of the curve fitting was observed for the third model under consideration ( N = 2 ,  Table 4. A more sophisticated version of the model with three Voce terms ( N = 3 ) and three backstresses ( M = 3 ) exhibited only a negligible improvement of the curve fitting. The approximation results are shown in Fig. 5 while the material parameter values are gathered in Table 5. The information about the residue (RSS-Residual Sum of Squares) is given in the figures.
To validate the methodology described above some analytical relationships were derived. The relations listed below hold for the Ch-R elastoplastic model and any uniaxial tension/compression process. They follow from the tensor equations which were given in Sect. 2.
According to Eqs. (3) and (17) for | − X |< k + R no plastic flow occurs and the stress can be determined from the linear relation which follows from Eq. (1), i.e.
If the material is actively yielding, Eq. (3) leads to the formula for the stress: with According to Eq. (8) the total backstress is given as: After integrating Eq. (9) and applying the initial conditions, it is found that: where p 0 is the initial plastic strain and X (i) 0 is the initial backstress value ( i = 1, 2, … , M ). In this particular case the initial values are interpreted as the accumulated values of plastic strain and backstress at the start of the current flow process.
It follows from Eq. (7) that the isotropic hardening is determined by N Voce terms, i.e.
with the effective plastic strain given as: The predictions of the analytical Eqs. (24-30) were used to verify the results obtained using the discretized version of the model, i.e., Eqs. (18)(19)(20)(21)(22). In Fig. 6, the stress values calculated for an exemplary cyclic tension/compression process are compared. It is seen that the predictions of the analytical equations and the numerical computations are in a prefect agreement. The validation test was performed for  Table 4.

Approximation of cyclic torsion-reverse torsion data
The general relations of the discretized model which are given in Sect. 3 were used to derive a set of equations describing the cyclic torsion-reverse torsion process. These equations are discussed below. They were utilized for calculating the theoretical shear stress values during the least squares optimization procedure.

Calculating theoretical stress values
In the case of torsion/reverse torsion process the stress tensor, its deviator and the backstress tensor have the following components: with X being the shear backstress component. For the computational step n + 1 the predictor shear stress is given as with Δ s being the shear strain increment. The plastic f low does not occur when the condition √ 3 | pr n+1 − X n |< k + R(ē p n ) is satisfied. In this case it is assumed that n+1 = pr n+1 . If the material is yielding the effective plastic strain increment Δē p needs to be calculated from the following nonlinear equation: where After solving Eq. (33) the shear stress value is updated, i.e.,

Identification of material parameters
The discretized model equations obtained for the torsion/reverse torsion process were used for computing the .
theoretical shear stress values during the material parameter identification procedure. The experimental data from the cyclic torsion-reverse torsion tests performed on copper were utilized. The total square error function took the form: with being the theoretical shear stress value calculated using Eqs. (32-35) and ̃ being the experimental stress measurement, whereas j is the number of the collocation point ( j = 1, 2, … , n ). Again, the approximation was performed in the Scilab software using the fminsearch function [22]. The material parameters were identified for the following versions of the elastoplastic model: The curve fitting results that were achieved for the simplest version of the Ch-R elastoplastic model ( N = 1 , M = 1 ) are shown in Fig. 7. The determined material parameters have been gathered in Table 2. Only slight improvement of the approximation quality was obtained in the case of the Yoshida model. Thus, the plot presenting the curve fitting is not included here, whereas the identified parameters of the model by Yoshida are listed in Table 3.
Adding an additional isotropic term and a backstress to the model ( N = 2 , M = 2 ) resulted in a significant improvement of the curve fitting quality, see Fig. 8. The identified material parameter values have been gathered in Table 4.
The model utilizing three isotropic terms ( N = 3 ) and three backstresses ( M = 3 ) achieved only negligible improvement of the approximation. The curve fitting results can be viewed in Fig. 9, whereas the determined parameter values are listed in Table 5.
All the approximations described above were performed using the same values of the Young's modulus, the Poisson's ratio and the initial yield stress as in the case of uniaxial tension/compression. The information about the residue (RSS) is provided in the figures showing the curve fitting.
For the purpose of verifying the results produced by the discretized version of the model given by Eqs. (32-35), the analytical equations describing the cyclic torsion/reverse torsion process were derived. They are attached below.
The plastic flow does not occur if the condition √ 3 | − X |< k + R is satisfied. In this case, the shear stress is calculated from a linear relation which follows from Eq.
where p s is the plastic shear strain. If the plastic flow occurs, the stress is calculated from the formula which can be derived from Eq. (3): The total shear backstress is given as a sum After integrating Eq. (9), the following relation for the backstress components is obtained: Fig. 7 Comparison of experimental results obtained for copper (cyclic torsion-reverse torsion) and theoretical predictions of Chaboche-Rousselier elastoplastic model with one isotropic hardening term ( N = 1 ) and one backstress variable ( M = 1) where p s0 is the initial plastic shear strain and X (i) 0 is the initial shear backstress value ( i = 1, 2, … , M ). In this particular case, the initial values are interpreted as the accumulated values of the plastic strain and the backstress at the start of the current flow process.
The isotropic hardening is governed by N Voce terms, i.e., with the effective plastic strain given as The predictions of the analytical Eqs. (38-43) were used to validate the results obtained using the discretized version of the model, i.e., Eqs. (32-35). In Fig. 10 the stress values calculated for an exemplary cyclic torsion/reverse torsion process are compared. It is seen that the predictions of the analytical equations and the numerical computations are in an excellent agreement. The validation test was performed for the material parameter values determined from the torsion/ reverse torsion process for N = 2 and M = 2 , cf Table 4.

Conclusions
In this work, a new method of determining the material parameters of mixed hardening elastoplasticity is presented. The method is applicable to the constitutive model proposed by Chaboche and Rousselier (Ch-R) [7,8] and its modifications. The presented method is highly automated and requires only basic processing of the experimental data. Furthermore, the method can be utilized for any loading histories.
The developed method was used for identifying the material parameters of copper. For that purpose, cyclic loading tests were performed on copper specimens that included uniaxial tension/compression and torsion/reverse torsion processes. The determined material parameter values of copper are gathered in Tables 1, 2, 3, 4, 5.
It is demonstrated that the Ch-R model utilizing two Voce terms to describe the isotropic plastic behavior and two backstress variables to simulate the kinematic hardening proved to provide a very good description of the elastoplastic response of copper (Figs. 4 and 8). Employing additional terms to capture the isotropic and kinematic hardening behaviors ( N = 3 , M = 3 ) results in a negligible improvement of the approximation quality (cf Figs. 5 and 9). Thus, it can be concluded that in the case of copper, the model with two Voce terms and two backstress variables is optimal and there is no reason in model's further extension within the discussed formulation of constitutive equation.
A certain shrinkage of the hysteresis loop was observed during the cyclic loading/unloading experiments that were performed on copper, e.g., Fig. 4. This phenomenon is interpreted as the isotropic softening effect. The negative parameter values that were identified for the Voce terms are responsible for capturing this effect (see Tables 2, 3

, 4, 5).
It is demonstrated that a model calibrated using the cyclic torsion data is able to properly describe the material response in tension/compression, cf Fig. 11a. On the other hand, it is shown that when the material parameters are determined using the tension/compression data the model is able to properly reproduce the material behavior in torsion (Fig. 11b). Thus, it can be concluded that for the considered group of constitutive equations the material parameter values can be identified using a single mechanical test (cyclic tension/compression or cyclic torsion). For the majority of the analyzed versions of the Ch-R elastoplastic model the material parameter values identified from the tension/ compression test are similar to those determined from the torsion test, cf Tables 2, 3, 4. However, for the model using three Voce terms to describe the isotropic hardening ( N = 3 ) and three backstress variables ( M = 3 ) some substantial differences between the identified parameter values can be observed. This fact can be explained by the growing complexity of the parameter value nonlinear optimization problem. This complexity is manifested in the case of the extended version of the model.
The correctness of proposed methodology has been validated using the derived analytical relations. An excellent agreement between the numerical and the analytical results was observed, cf Figs. 6 and 10.

Conflicts of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 11
Comparison of experimental results obtained for copper and theoretical predictions of Chaboche-Rousselier elastoplastic model with two isotropic hardening terms ( N = 2 ) and two backstress variables ( M = 2 ): a uniaxial test data and model predictions for mate-rial parameters determined from torsion test, b torsion test data and model predictions for material parameters determined from uniaxial test