Cyclic Elasto-Plastic Behaviour of 51CrV4 Steel for Railway Parabolic Leaf Spring Design

1. Introduction

Over the years, engineers and researchers have developed more efficient geometries and processes to maximize the strength of suspension elements. Although leaf springs have been elements used for many years, these components have evolved both in terms of their resistance and geometries. The evolution of this component is also dependent on the type of sectors where these elements are inserted, as for example, road vehicles [1,2,3,4,5,6,7,8,9,10] and rail vehicles [11,12,13,14,15]. In the case of the railway freight sector, parabolic leaf springs (illustrated in Figure 1)) that are standardized by the International Union of Railways (UIC) [16] have gained higher representation due to their optimized geometry, which allows an increase in the specific energy storage capacity and reduction of the friction effects, which consequently increases the safety and quality of the wagon motion [17,18].

This two-axle freight vehicle suspension component is produced from high-strength steel, designated as 51CrV4 [19,20,21,22], which might be applicable in various industry machinery components in addition to road and rail suspension vehicles [23,24], such as steering knuckles, connecting rods, spindles, pumps, and gear shafts, pinions, forged crankshafts, and gears [21,25,26].

Despite the good mechanical resistance characteristics presented in several investigations [24,27,28,29], due the presence of geometric notches [18,30,31,32,33] and high surface roughness of its leaves (average up to $R_a$ = 7 $\mu$m [34] and maximum roughness, $R_z$ = 40 $\mu$ m [35]), local stress levels can reach or exceed the yield limit of the material. Operating under these cyclic loading conditions, the cyclic elasto-plastic material behaviour model must be considered.

Over the years, different constitutive models have been developed to represent cyclic plastic behaviour. Some examples are the linear kinematic hardening model proposed by Prager [36], the multiple-yield surface model [37], two-surface plasticity model, yield and boundary surface model, [38,39], and others such as Malinin and Khadjinsky [40], Armstrong and Frederick [41], Chaboche [42], and Jiang [43]. In fact, the development of the cyclic plasticity models helped in the comprehension of the different mechanical responses associated with the plastic domain. Transient effects associated with Bauschinger’s effect, ratcheting, hardening or softening straining, as well as stress relaxation, are some of the examples observed and investigated more frequently in cyclic plasticity [42,44,45,46,47,48,49,50,51,52,53,54,55,56,57].

Some of these cyclic plasticity models have been used in numerical models to perform fatigue prediction analysis in cyclic elasto-plastic conditions for both smooth and notched components [58,59,60,61,62,63]. Even in components subjected to cyclic loadings below the yielding strength, due to the surface roughness, the localised plastic deformation may occur [64,65,66]. Under these conditions, cyclic plasticity should be considered for fatigue crack initiation analysis. Moreover, in fatigue crack propagation analysis, there is localised plasticity at the crack tip. Thus, plasticity models should also be considered in determining the cyclic material behaviour at crack tip [67,68,69,70].

In this investigation the Chaboche model, widely used for designing structural components under cyclic loading, is considered. The Chaboche model considering the effect of combined kinematic and isotropic hardening is calibrated using the experimental data on 51CrV4 martensitic steel, obtained in strain-controlled fatigue tests. This combination of kinematic and isotropic hardening models has been widely used in different works, resulting in good prediction results compared with experimental data [45,71,72,73]. The parameters of the Chaboche kinematic hardening model are calibrated considering the stabilised cyclic elastoplastic Ramberg-Osgood curve for 90% of the total life of the specimens. Regarding the isotropic hardening model parameters, these are calibrated considering the material softening amount curve ([74]) and the curve of variation in the size of the normalized elastic limit surface throughout the test. In the parameter calibration procedure, the Hager-Zhang’s conjugate gradient unconstrained optimization method was considered. To validate the estimated combined hardening parameters, the estimated parameters are used in the numerical computation using an elasto-plastic finite element approach. The outcomes by numerical model are subsequently compared with experimental data and the comparison errors are determined.

2. Cyclic Elasto-Plasticity for Steels

2.1. Elasto-Plasticity Modelling

Cyclic plasticity, a crucial factor in fatigue analysis, is examined for components subjected to loads exceeding the material’s yield strength. The Chaboche hardening model [42] is used to characterize this cyclic behaviour. This model is based on the yield von Mises plasticity criterion, as outlined below

$$F(\sigma -\alpha ,{\sigma }_y)=\sqrt{({\sigma }^{\prime }-{\alpha }^{\prime }):({\sigma }^{\prime }-{\alpha }^{\prime })}-({\sigma }_y+R)=0,$$

(1)

with ${\sigma }_y$ and R denoting, respectively, the initial size of the elasticity limit surface and the change in the size of the elasticity limit surface. If $R>0$, the elasticity limit surface expands, otherwise, F contracts. Concerning the stress tensors, ${\sigma }^{\prime }$ and ${\alpha }^{\prime }$ denoting respectively the deviatoric part of the stress tensor, $\sigma$, and the back stress tensor, $\alpha$.

The direction of plastic flow is determined by the associative flow rule, which follows the principle of normality, such as

$${\dot{\epsilon }}^p=\dot{\lambda }{\displaystyle \frac{\partial F}{\partial \sigma }},$$

(2)

where ${\dot{\epsilon }}^p$ denotes the time rate of plastic strain tensor. The consistency condition, $\sum \left[\partial f/\partial \theta \right]\dot{\theta }=0$, serves as the condition for updating the stress and strain-state variables, with $\theta$ assuming $\sigma$, $\alpha$, or R. The plastic multiplier in Equation (2), $\dot{\lambda }$ is determined from the consistency condition, which ensures the Kuhn–Tucker conditions.

2.2. Kinematic Hardening Model

In several investigations, the cyclic hardening behaviour, mainly in steels, has been modelled by the Chaboche hardening model. Chaboche kinematic hardening was developed considering the von Mises plasticity criterion (Equation (1)). As advantage, Chaboche model can predict the transient effects such as cyclic hardening or softening, ratcheting, mean stress relaxation, etc. Chaboche kinematic hardening model is an improved model from Amstrong and Frederick, AF model, which itself is an improved model from the linear hardening rule [36], developed by Prager [36], by introducing the concept of dynamic recovery [41] . Thus, the back stress tensor, $\alpha$, is given by a sum of i-th component of back stress tensor, ${\alpha }_i$, such that [42]

$$\alpha =\sum _{i=1}^N{\alpha }_i,$$

(3)

and the variation of each individual back stress tensor, ${\dot{\alpha }}_i$, is

$${\dot{\alpha }}_i={\displaystyle \frac{2}{3}}{C}_i\dot{{\epsilon }^p}-{\gamma }_i{\alpha }_i\dot{p}.$$

(4)

In equations (3) and (4), $C_i$ and ${\gamma }_i$ are fitting parameters and $\dot{\mathbf{p}}$ denotes the effective plastic strain rate, which is given by

$$\dot{\mathbf{p}}=\sqrt{{\displaystyle \frac{2}{3}}{\dot{\epsilon }}^p:{\dot{\epsilon }}^p}.$$

(5)

The Chaboche kinematic hardening model permits overcome some constraints pointed out by Prager and AF models, such as the difficulty of a accurate behaviour representation for large deformation regimes, and the transition from elastic to plastic regimes [42]. Several investigations have proved that considering only three components can reach good predictions ratcheting behaviour [42,55].

2.3. Isotropic Hardening Model

Although the Chaboche hardening model is modeled for material behaviour with kinematic hardening, it can be considered with isotropic hardening. The isotropic hardening that represents the change in the size of the elasticity limit surface, $\dot{R}$, is given by

$$\dot{R}=b_{\infty }\left(R_{\infty }-R\right)\dot{p},$$

(6)

where $b_{\infty }$ denotes the rate of isotropic hardening and $R_{\infty }$ denotes the saturation of isotropic hardening. Both $b_{\infty }$ and $R_{\infty }$ are parameters depending on the material and temperature. Lower values of $b_{\infty }$ indicate a slower saturation of the model, whereas for materials exhibiting a high value of $b_{\infty }$, its yield surface will to evolve very fast cycle-by-cycle. Integrating the Equation (6) results that

$$R=R_{\infty }\left[1-exp\left(-b_{\infty }p\right)\right].$$

(7)

Equation (7) provides the size of the elasticity limit surface, which may be considered for both in monotonic and cyclic loading conditions.

2.4. Empirical Estimation of Hardening Variables

The estimation of kinematic and isotropic hardening variables is normally carried out considering test data under strain-controlled cyclic loading conditions. Under these conditions, the kinematic hardening is represented by equation (8) as [42]

$${\sigma }_a={\sigma }_y^{\prime }+{\displaystyle \frac{\Delta \alpha }{2}}=\left({\sigma }_y+R_y^{\prime }\right)+\sum _{i=1}^{n_{\alpha }}{\displaystyle \frac{C_i}{{\gamma }_i}}tanh\left({\gamma }_i{\epsilon }_a^p\right),$$

(8)

where the parameters $C_i$ and ${\gamma }_i$ are determined by fitting the Equation (8), and $R_y^{\prime }$ denotes the variation of the elasticity limit surface for the stabilized-state hysteresis loop. Equation (8) allows taking into account both the ratcheting effect and the stress saturation effect [42,75]. According to some research, the limit stress stabilised cyclic curve is reached for strain amplitudes of 2% and 3% in high-strength steels, which means that after 2% and 3%, the increase of stress amplitude with the strain amplitude is quite low [76,77,78]. This slight increase of hardening stress may be considered summing to the equation (8) a linear component of stress, given by Prager’s model, which results in

$${\sigma }_a=\left({\sigma }_y+R_y^{\prime }\right)+\sum _{i=1}^{n_{\alpha }-1}{\displaystyle \frac{C_i}{{\gamma }_i}}tanh\left({\gamma }_i{\epsilon }_a^p\right)+{\displaystyle \frac{2}{3}}{C}_i{\epsilon }_a^p.$$

(9)

The fitting of the parameters $C_i$ and ${\gamma }_i$ can be carried out considering the hysteresis loop in stabilized conditions for the largest strain amplitude analysed, or by the elasto-plastic cyclic curve. According to the literature [42,74,79,80], the Ramberg-Osgood, RO, model is one of the most suitable model for fitting the Chaboche’s parameters in steel materials. The RO model is a stable elastic-plastic hardening that relates the total strain amplitude, ${\epsilon }_a^t$, to the stress amplitude, ${\sigma }_a$, through a power law, of the form [81], such that

$${\epsilon }_a^t={\epsilon }_a^e+{\epsilon }_a^p={\displaystyle \frac{{\sigma }_a}{E}}+{\left({\displaystyle \frac{{\sigma }_a}{K^{\prime }}}\right)}^{{\displaystyle \frac{1}{n^{\prime }}}},$$

(10)

where E is the Young’s modulus, $K^{\prime }$ is the cyclic hardening strength coefficient, and $n^{\prime }$ is the cyclic strain hardening exponent. Ramberg–Osgood’s model (Equation (10)) is given by decomposing the total strain amplitude, ${\epsilon }_a^t$, into its elastic, ${\epsilon }_a^e$, and plastic, ${\epsilon }_a^p$, parts. In spite of being an ancient model, several researchers [24,74,82,83,84,85,86,87,88] recently used this equation model (10) to characterize the elasto-plastic behavior of a material under cyclic loadings. , Also, data obtained from strain-controlled cyclic testing, the parameters of the evolution of the elasticity limit surface (equation (7) can be obtained. Thus, the elastic limit surface (equation (7) can be represented in terms of the normalized maximum stress, $R^{*}$ as a function of the accumulated plastic strain, such that

$$R^{*}={\displaystyle \frac{R}{R_{\infty }}}={\displaystyle \frac{{\sigma }_{max}\left(N_i\right)-{\sigma }_{max}\left(N_{i=1}\right)}{{\sigma }_{max}\left(N_{i=s}\right)-{\sigma }_{max}\left(N_{i=1}\right)}},$$

(11)

with ${\sigma }_{max}\left(N_1\right)$, ${\sigma }_{max}\left(N_{i=s}\right)$, and ${\sigma }_{max}\left(N_i\right)$ denoting the maximum stress corresponding to the first cycle, the stabilized cycle, and the i-th cycle, respectively. $R_{\infty }$ denotes the difference between the value of yielding stress in the stabilized cycle and the value of yielding stress in the first cycle. $R_{\infty }$ relies on the imposed strain amplitude. Increasing the total strain amplitude, the amount of hardening/softening, $R_{\infty }$, increases with the maximum value, $R_{\infty ,max}$. In several investigation [76,77,78], it was found that for medium and high strength steels, the variation in the elasticity limit surface remain unchanged for around 2 or 3 % of applied strain. $R_{\infty }$ varies with the applied strain amplitude, ${\epsilon }_a$, through a decreasing power law, written as [89]

$$R_{\infty }=R_{\infty ,max}\left[1-exp\left(k_0\left({\epsilon }_a-{\epsilon }_0\right)\right)\right],$$

(12)

where parameters $k_0$ and ${\epsilon }_0$ are estimators for the amount of cyclic softening. This approach can be considered because increasing the total strain amplitude, the amount of hardening/softening, $R_{\infty }$, tends to increase up to the maximum value, $R_{\infty ,max}$. $R_{\infty ,max}$ can be estimated for applied strain amplitude values around 2% or 3% for medium and high strength steels [76,77,78].

Concerning the hardening parameter that governs the saturation rate, $b_{\infty }$, this parameter may also be computed explicitly as

$$b_{\infty }=-{\displaystyle \frac{ln\left[1-R^{*}\right]}{2N_f\Delta {\epsilon }^p}},$$

(13)

according to the methods described in [42,90,91].

3. Material and Analysis Procedure

The calibration of the cyclic behaviour of the steel under investigation, 51CrV4, is carried out considering two approaches, experimental and computational. In the experimental approach, data collected from cyclic tests are used to estimate kinematic and isotropic hardening parameters, using statistical fitting models. The estimated parameters are considered in the computational models in order to compare and validate the Chaboche model for 51CrV4 spring steel applied to parabolic leaf springs of railway freight wagons.

Table 1. Standard chemical composition of 51CrV4 steel grade in % wt .

Material	C	Si	Mn	Cr	V	S	Pb	Fe
DIN 51CrV4 (1.815)	0.47-0.55	≤0.40	0.70-1.10	0.90-1.20	≤0.10-0.25	≤0.025	≤0.025	96.45-97.38

Table 1. Standard chemical composition of 51CrV4 steel grade in % wt .

Material	C	Si	Mn	Cr	V	S	Pb	Fe
DIN 51CrV4 (1.815)	0.47-0.55	≤0.40	0.70-1.10	0.90-1.20	≤0.10-0.25	≤0.025	≤0.025	96.45-97.38

3.1. Chemical Composition and Microstructure

51CrV4 steel is a chromium-vanadium alloyed steel with an average carbon content of 0.50%, approximately. When used in parabolic leaf springs, 51CrV4 steel follows the heat treatments described in ISO 683-14 standard [22], which consists of quenching at 850 ^oC in an oil bath and then tempering at 450 ^o W . As a result, 51CrV4 presents a microstructure of quenched martensite with retained austenite (white phases) as shown in Figure 2[92].

The heat treatment provides the 51CrV4 steel its high mechanical strength. According to the material characterization tests carried out in [24] under monotonic loading conditions, following the ISO 6892-1 [93] standard, the spring steel exhibits a yield strength, ${\sigma }_y$ = 1271.48 MPa, an ultimate tensile strength, ${\sigma }_{uts}$ = 1438.5 MPa, a (conventional) ductility at fracture, ${\epsilon }_f$ = 7.53 %, and a Vickers hardness of 447 HV. Table 2 presents the summary of the 51CrV4 mechanical properties, with E denoting the Young’s modulus and $R_{A_f}$, the area of reduction at fracture.

3.2. Experimental Procedure and Statistical Techniques

Kinematic and isotropic mechanical properties are to be determined for chromium-vanadium alloyed steel. The cycle tests carried out in [24] are considered to determine the cyclic hardening properties and subsequent verification in the computational model by analysing the hysteresis loops. The sample set consists of seven flat-shaped specimens with a net cross-section, $A_0$, of 38.98 ± 0.28 mm² whose width, $w_o$, is $6.04\pm 0.04$ mm and thickness, $h_o$, is $6.46\pm 0.01$ mm (see Figure 3 - B) (see Table 3).

The specimens were manufactured and tested under constant strain-controlled amplitude loading conditions according to ASTM E606 standard [94] for a strain ratio, $R_{\epsilon }=0.0$ and an average strain rate d$\epsilon$/dt of 0.8 %/s. The applied strain amplitudes were 0.375 ($\times 2$), 0.50 ($\times 2$), 0.750, 1.00, and 1.25 %. The cyclic tests were conducted in an INSTRON 8801 servo-hydraulic machine with a load cell of 100 kN, with a dynamic INSTRON 2620-202 clip gauge with a working range of ± 2.5 mm (see Figure 4).

The kinematic and isotropic hardening properties are estimated using statistical techniques based on the least-squares method. In the case of the problem can be described by a linear function, the ordinary least-squares (OLS) is considered, whereas if the problem is non-linear, an unconstrained numerical optimization scheme is assumed.

For a linear response function, then

$$\mathbf{y}={\beta }_0+{\beta }_1\mathbf{x},$$

(14)

with $\mathbf{x}$ and $\mathbf{y}$ denoting the design vector and respective response, respectively [95]. The minimization problem results in the equation system in (15), with

$$\widehat{\beta }=arg\underset{\widehat{\beta }\in \Omega }{min}{\left|\right|\mathbf{y}-{\widehat{\mathbf{y}}}_i\left|\right|}^2={\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }\mathbf{y},$$

(15)

where $\mathbf{X}$ is the design matrix, and ${\mathbf{X}}^{\prime }$ is its transposition. $\widehat{\beta }$ is the estimator vector containing estimators, which in this case is ${\beta }_0$ and ${\beta }_1$.

On the other hand, for a non-linear response function, the objective function

$$\widehat{\beta }=arg\underset{\widehat{\beta }\in \Omega }{min}{\left|\right|\mathbf{y}-{\widehat{\mathbf{y}}}_i\left|\right|}^2=\underset{\widehat{\beta }\in \Omega }{min}\left[{\displaystyle \frac{1}{N_{\mathbf{x}}}}\sum _{i=1}^{N_{\mathbf{x}}}{\left({\displaystyle \frac{{\mathbf{y}}_i-{\widehat{\mathbf{y}}}_i}{{\widehat{\mathbf{y}}}_i}}\right)}^2\right],$$

(16)

with $N_{\mathbf{x}}$ denoting the size of vector $\mathbf{x}$ is assumed. The minimization numerical problem is solved using the iterative method of Broyden–Fletcher–Goldfarb–Shanno (BFGS), with the recurrence equation given by [96]:

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k^s{\mathbf{s}}_k,$$

(17)

where in each iteration, k, the search direction, ${\mathbf{s}}_k$, is obtained using Hager-Zhang’s conjugate gradient method [97], for a maximum convergence error of 1.00 × 10⁻⁵. The gradient of the objective function for x is computed using forward-mode automatic differentiation [98] and the step size, ${\alpha }_k^s$, is always positive and assumed to be 1.0. With respect to the initial guess for a hessian matrix ${\mathbf{H}}_0$, this one is the identity matrix.

3.3. Computational Approach

A computational approach using finite element analysis (FEA) in ANSYS software was adopted to verify the cyclic material properties derived from the statistical analysis of the experimental data. A pure-displacement finite element formulation was considered for a cyclic small elasto-plastic deformation analysis of the material. The system of differential equations was solved incrementally by applying a displacement increment, $\Delta \delta$, according to the Newton-Raphson procedure (with the L2-Norm), and using the direct sparse method with LDLT factorization for the matrix equation (18),

$${\mathbf{K}}_t\Delta \mathbf{u}={\mathbf{F}}_{ext}-{\mathbf{F}}_{int},$$

(18)

with ${\mathbf{F}}_{ext}$ denoting the external forces applied at nodes, ${\mathbf{F}}_{int}$ representing the restoring loading vector, and ${\mathbf{K}}_t$ denoting the symmetric and positive-definite global tangent matrix.

A single quadratic 20-node iso-parametric finite element with dimensions $1\times 1\times 1$ following a full integration scheme of Gauss was assumed. The von Mises plasticity criterion combined with the Chaboche (kinematic plus isotropic) hardening model [42] was considered in the integration of the plasticity law, via the return mapping algorithm, according to the Euler backward implicit scheme [99,100,101]. The finite element simulates the material behaviour under a uniformly distributed strain condition, by considering isostatic boundary conditions.

4. Results and Discussion

4.1. Cyclic Behaviour

Cyclic tests were performed until failure as illustrated in Figure 5. In research carried out in [24], it was found that 51CrV4 steel exhibits cyclic softening behaviour. Figure 5 shows the hysteresis loops for the 1st cycle and the cycles corresponding to 50 and 90 % of the total life. The consideration of the hysteresis loop for 50% of the total life is widely used for the development of fatigue resistance curves, however the 51CrV4 steel revealed in the analysis in [24] that for high strain amplitude conditions (higher strain amplitudes than 0.5 %), the hysteresis loops do not stabilise completely, which is also proven by Figure 5, with a continuous cyclic softening and mean stress relaxation, however the differences and variations are quite low. Thus, in the analysis of the elasto-plastic transient given by the isotropic and kinematic cyclic laws for a life of 90% of the total life is considered. Assuming 90 % of the total life, it is still ensured that the hysteresis loops do not consider the final loops containing crack initiation and propagation phenomena.

4.2. Determination of the Kinematic Hardening Properties

When determining the kinematic hardening parameters, hysteresis loops for a life of 90% are considered. The hysteresis loops for a 50% of the total life are also used considered, but only for comparison. Thus, considering the RO cyclic curve (equation (10)) for 90% of total life, cyclic RO parameters of $K^{\prime }=1576.2$ MPa, and $n^{\prime }=0.0918$ with a coefficient of determination, $R^2=0.8113$ are obtained. The kinematic hardening properties are determined by fitting the curve of equation (9) to the cyclic RO curve (equation (10)), using the unconstrained numerical optimization scheme previously described, where in equation (16), such that

$$\widehat{\mathbf{y}}=\sum _{i=1}^{n_{\alpha }-1}{\displaystyle \frac{C_i}{{\gamma }_i}}tanh\left({\gamma }_i{\epsilon }_a^p\right)+{\displaystyle \frac{2}{3}}{C}_i{\epsilon }_a^p$$

and

$$\mathbf{y}=K^{\prime }{\left({\epsilon }_a^p\right)}^{n^{\prime }}.$$

with number of back stress, $n_{\alpha }=1,2,3,4$. The selected model was the model that provided the highest minimization of the objective function, hence the smallest $RRMSE$, the smallest number of back stress, and that ensured the condition, $C_i$ and ${\gamma }_i>C_{i+1}$ and ${\gamma }_{i+1}$ with $i=1,..,n_{\alpha }$. As a result of this condition, the selected Chaboche’s kinematic model is constituted by 4 back stresses (3 non-linear and 1 linear).

Since the optimal solution is dependent on the input parameters, the process to obtain the best initial design vector, $\mathbf{x}=\left\{C_i,{\gamma }_i\right\}$, with $i=1,...,n_{\alpha }$, is made from the generation of a random sample. The generation of the random values is made using a pseudo-generator of random values designed as Mersenne Twister that is a twisted GFSR random number generator of 32 bits MT19937 [101,102]. Side bounds are introduced in the model to constrain the regression parameters $C_i$ and ${\gamma }_i$. The lower and upper side bounds are assumed to be respectively for both parameters, ${10}^{-7}$ and ${10}^7$. Parameters are fitted until the value of fracture strain, ${\epsilon }_f=7.53\%$. Notice that the fourth back stress is linear in order to take into account the effect of ratcheting and mean stress relaxation effect [45].

Considering the Chaboche parameter fitting method for 4-back stress kinematic hardening models (3 nonlinear and 1 linear), the optimal parameters are $C_1=81701.2$ MPa, ${\gamma }_1=610.2$, $C_2=22035.4$ MPa, ${\gamma }_2=159.14$, $C_3=6332.2$ MPa, ${\gamma }_3=39.93$, and $C_4=2673.4$ MPa, with a relative root mean square, $RRMSE$, from 0.2853 %. The value of each fitting parameter, $C_i$ and ${\gamma }_i$ is presented in Table 4.

Figure 6 depicts the decomposition of total the back stress tensor, ${\alpha }_t$, in the model with 3 non-linear, ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$ and 1 linear, ${\alpha }_4$ components. From Figure 6 - Right, it is notable that there is a good agreement between the variation of back stress components with the the assumptions of Chaboche’s model. ${\alpha }_1$ and ${\alpha }_2$ have higher impact in the initial hardening until 2 % of the plastic strain amplitude. On the other hand, ${\alpha }_3$ contributes gradually to increase the kinematic hardening from the beginning of the plastic deformation process. Lastly, ${\alpha }_4$ increases linearly the total back stress tensor. In addition, Figure 6 - Left compares the decomposition of total stress, $\sigma$, into the total back stress, ${\alpha }_t$, and the size of elasticity limit surface.

In case hysteresis loops of 50 % of the total life were chosen, the RO cyclic curve described with $K^{\prime }=1514.3$ MPa (more 4.09 %), $n^{\prime }=0.0790$ (more 16.20 %), and ${\sigma }_y^{\prime }$ of 731.5 MPa (less 7.49 %) results in $C_1=80097.1$ MPa, ${\gamma }_1=648.4$, $C_2=20486.2$ MPa, ${\gamma }_2=164.9$, $C_3=5650.1$ MPa, ${\gamma }_3=40.64$, and $C_4=2673.4$ MPa (presented in Table 4), whose highest differences occur for constants $C_3$ and $C_4$. Note that in the case where the RO cyclic model parameters were fitted by estimation models using monotonic properties, for example the Modified Universal Slopes method, MUSM, (with coefficients ${\sigma }_f^{\prime }=2057.36$ MPa, ${\epsilon }_f^{\prime }=0.2339$ exponents $b_f^{\prime }$ = -0.09 and $c_f^{\prime }$ = -0.56 [103]), results in $C_1=82787.4$ MPa, ${\gamma }_1=580.07$, $C_2=22260.6$ MPa, ${\gamma }_2=136.5$, $C_3=5835.24$ MPa, ${\gamma }_3=31.29$, and $C_4=2305.5$ MPa, with a $RRMSE$ of 0.3167%, and a cyclic yield strength ${\sigma }_y^{\prime }=579.1$ MPa (which results in a variation in the size of the elastic limit surface, $R_y^{\prime }$ = -510.1 MPa). Comparing both Chaboche kinematic data, it is verified that the constants C are in the same order of magnitude, except the constants $\gamma$. These constants obtained from the experimental data at $R_{\epsilon }=0$ tend to be lower than $R_{\epsilon }=-1$ (MUSM method). The highest reduction occurred in ${\gamma }_1$, which is responsible for the first hardening zone, which along the way is the zone most affected by the parameter $n^{\prime }$.

4.3. Determination of the Isotropic Hardening Properties

The isotropic hardening law properties, $R_{\infty }$ and $b_{\infty }$, in equation (6) are determined separately. Initially the parameter $R_{\infty }$ is determined and subsequently $b_{\infty }$ is determined, both using the unconstrained optimization scheme.

Considering the determination of the isotropic hardening parameters following the unconstrained optimization scheme for the variation of $R_{\infty }$, (equation (12)), then, in accordance with equation (16), it is assumed that

$${\widehat{\mathbf{y}}}_i=R_{\infty ,max}\left(1-exp\left[k_0\left({\epsilon }_{a,i}-{\epsilon }_0\right)\right]\right)$$

and

$${\mathbf{y}}_i=R_{\infty ,i}.$$

In the numerical procedure, the lower and upper bounds of regressors, $R_{\infty ,max}$, $k_0$, and ${\epsilon }_0$ are found. After finding the side bounds, a pool of permissible solution of regressors is created in order to be considered as the initial design vector, $\mathbf{x}=\left\{R_{\infty ,max},k_0,{\epsilon }_0\right\}$. From a pool of optimal solutions, the best optimal solution is considered. It is assumed as the best optimal solution, the solution that has the lowest value of the objective function. According the algorithm, the optimal solution of $\mathbf{x}$ is $R_{\infty ,max}=-232.2$ MPa , $k_0$ = - 380.55, and ${\epsilon }_0=0.32\%$, with a relative root mean square error, $RRMSE=14.83\%$.

Figure 7 depicts the evolution of the amount of cyclic softening as a function of the applied strain amplitude. According to this variation, the existence of a cyclic softening is considered after a threshold strain value, ${\epsilon }_0$. The threshold strain, ${\epsilon }_0$ corresponds to the strain for a stress of 648 MPa, which is approximately equal to the cyclic yield strength, ${\sigma }_{01}$ = 676.7 MPa, but 4.24 % inferior. From Figure 7, it is verified that $R_{\infty ,max}$ is reached after 2.0 % of strain amplitude, which is in accordance with [76,77,78].

Regarding the isotropic saturation rate parameter, $b_{\infty }$, its value is determined according to the method in [42,90,91], such that the response curves are assumed to be

$${\widehat{\mathbf{y}}}_i=1+exp(-{\widehat{b}}_{\infty ,i}{\mathbf{p}}_i)$$

and

$${\mathbf{y}}_i=R_i/R_{\infty }.$$

Figure 8 depicts the normalized variation of the yield surface size evolution, $R/R_{\infty }$, throughout the cyclic testing until 90 % of the total life to failure. Both the mean curve and standard deviation bands concerning the mean value of $b_{\infty }$ are also presented. From the comparison between experimental data and numerical outcome, it is notable that cyclic tests performed at ${\epsilon }_a=0.375\%$ exhibited a no-coherent evolution, due to its large differences in rate size evolution. In contrast, cyclic tests carried out at higher strain amplitudes than $0.375\%$ exhibited a better agreement with respect to the numerical approach.

Considering the 90% of total-life criterion, most cyclic tests appear not to have an actual stabilization plateau as suggested by numerical approach [42]. However, one of the cyclic tests, at ${\epsilon }_a=0.375\%$ appears to have a stabilised state. With respect to the quality of the fitting curve, it is noted that in the initial cyclic softening phase, from ${10}^{-2}$ until 1, the regression curve tends to exhibit a lower saturation rate whereas in the middle phase, the softening evolution follows the same trend than the experimental results. This behaviour has been also found in some materials [90,104]. The average value found for $b_{\infty }$ was 0.3705 $\pm 0.1628$ with a $RRMSE$ of $42.00\pm 13.22\%$. According to [90], $b_{\infty }$ is usually assumed to be between 50 and 0.5 considering the half-life criterion. In spite of $b_{\infty }$ being slightly lower than 0.5, $b_{\infty }$ represents the average value obtained from different tests for a full cyclic test. The high value of $RRMSE$ is related to the shape of the softening evolution curve. In order to reduce the $RRMSE$ value, other isotropic hardening models contemplating an improvement in the stabilization rate and the shape of the softening/hardening evolution curve may be considered [104,105].

4.4. Cyclic Elasto-Plastic Response Using a Computational Approach

After determining the best fitting for the kinematic and isotropic hardening parameters, the numerical model following the finite element method approach as described previously, with the respective elasto-plastic mechanical properties presented in Table 5 is considered and compared with the experimental results to validate the cyclic plasticity model of chromium-vanadium steel 51CrV4.

For each strain amplitude condition, enough number of cycles were applied such that the stabilised behaviour will be possible to observe. The adopted strain ranges were only 1.00, 1.50, 2.00, and 2.50 %. For strain ranges inferior than 1.00 % it is observed a elastic hysteresis loops after the first cycles, and therefore the strain level of 0.75% was not considered for the analyses.

In a first phase, the hysteresis loops of the material numerical model were evaluated by comparing the first cycle and the cycle associated with the loop corresponding to 90% of the total life. The higher strain ranges were considered for the first analysis. Figure 9 depicts the comparison between experimental data in the first and stable hysteresis loops with the numerical outcomes in terms of stress and applied total strain for the presented strain amplitudes.

From the analysis of the monotonic part, it is noted the existence of a lag in the initial hardening between experimental and numerical outcomes. This lag is associated with the definition of the isotropic hardening law, such that the linear hardening coefficient (as described in the Voce’s law) was assumed to be null and also due to the material plastic properties variability. On the other hand, the maximum stress for each applied loading level tends to be very close between the models as illustrated in Figure 9, Figure 9, and Figure 9. In the case of Figure 9, it is verified that the yield strength is much higher than the others, generating a higher difference associated with the hardening material. However, this difference tends to remain constant since the experimental and numerical curves tend to be collinear.

Comparing the hysteresis loops associated with the first cycle, it is visible that the plasticity model cannot describe correctly the lower branch of the experimental model. This plasticity model imprecision is verified in both the three highest applied strain levels. However, comparing the experimental data associated with the first and steady-stable loops is notable that the evolution of the minimum stress is much less than the maximum tensile stress. This behaviour is typical for material exhibiting a non-Masing behaviour [58]. The variation of maximum stress in each cycle for data under analysis is presented in Figure 10. This figure clearly shows the existence of a non-symmetry associated with the maximum and minimum stress in each cycle of the softening process. In fact, the variation of the minimum stress of the hysteresis loop seems to be only associated with the kinematic hardening component for any applied strain range. This minimum stress tends to remain unchanged throughout the life of the specimen regardless of the applied strain range.

Concerning the effect of the isotropic saturation rate parameter, $b_{\infty }$, it is shown that the plasticity model with an average value of $b_{\infty }$ equal to 0.37 cannot describe correctly the first phase of the material softening for the higher strain ranges. However, for a strain range of 0.5 %, there is a good correspondence between the experimental and numerical behaviour for the maximum stress in each cycle. Notice that $b_{\infty }$ is an average value assumed constant for several strain amplitudes and throughout the cyclic life. Therefore, despite this value not corresponding to the reality of the material model, $b_{\infty }$ = 0.37 appears to be conservative for most of the whole life as illustrated in Figure 10.

Now, analysing the hysteresis loops in the steady-state conditions, one verifies that the maximum and minimum stress of the hysteresis loops are very similar compared with experimental data. In general, for higher applied strain ranges, as presented in Figure 9 and Figure 9, the maximum stress tends to be around 6.0 % superior, whereas, for the testing condition in Figure 9, the maximum stress is slightly inferior, around -4.0 %, in face to the experimental data (see Figure 11). Concerning the minimum stress, the offset of the numerical models is around -10.0 % as seen in Figure 11. The maximum deviation occurs between 5 and 30 % of the total-life and they can achieve 16 %. For specimens tested with a strain amplitude of 0.5 %, a different error analysis needs to be made. In fact, analysing Figure 9 it is verified that from the first cycle until the stabilised cycle, the strain ratio of $R_{\epsilon }=0.0$ is never achieved. Over this assumption, the error associated with the minimum value of the stress at the stabilised cycle is normal being above 100 % as seen in Figure 10. Concerning the hysteresis loop area of $\sigma -\epsilon$, one verifies that the area obtained from numerical outcomes tends to be conservative.

5. Conclusions

The modelling of the elasto-plastic behaviour in cyclic conditions of the chromium-vanadium alloyed steel (DIN 51CrV4) was obtained in this paper. The experimental data concerning the cyclic tests performed previously in [24] were used for determining the parameters of the isotropic and kinematic hardening plasticity model. Chaboche’s kinematic hardening model along with isotropic hardening model were considered for modelling the cyclic behaviour. The hardening model parameters were obtained from statistical techniques based on the least-squares method. These parameters were considered in the finite element model whose outcomes were posteriorly compared to the experimental results.

Initially, the Chaboche’s kinematic hardening properties for stabilised cyclic conditions assuming 90 % of the total life criterion were determined. The Ramberg-Osgood’s elasto-plastic parameters were used for fitting the kinematic coefficients $C_i$ and ${\gamma }_i$. The best fit was found for three non-linear and one linear back stress, with values of $C_1=81701.2$ MPa, ${\gamma }_1=610.2$, $C_2=22035.4$ MPa, ${\gamma }_2=159.14$, $C_3=6332.2$ MPa, ${\gamma }_3=39.93$, and $C_4=2673.4$ with a $RRMSE$ of 0.2853 %. These values were posteriorly compared with the ones obtained from stabilised cyclic conditions for a 50 % of the total life since this criterion is usually used for determining the fatigue life in strain-controlled life. Comparing both criteria, one verified that both criteria may be used without existing large differences in the estimators $C_i$ and ${\gamma }_i$, although the cyclic strength tends to decrease. Additionally, Chaboche kinematic hardening properties were estimated from estimation models based on monotonic properties, the modified Universal Slopes (MUSM), which does no affected significantly the kinematic coefficients if the cyclic yield stress is defined correctly.

Concerning the isotropic hardening parameters, the evolution of the yield surface, R, was determined for several strain ranges. A relationship for estimation of $R_{\infty }$ as a function of the applied strain amplitude is assumed. According to the relative root mean square error, a regression model with $RRMSE$ less than 15 % was found. In addition, a limit value for $R_{\infty }$ was found ($R_{\infty ,max}$) from the regression model. This value, $R_{\infty ,max}$, is independent of the applied strain range and is equal to -232.2 MPa. Regarding the isotropic saturation rate, $b_{\infty }$, from the optimization methods was found an average value of 0.3705 with a $RRMSE$ = 42.00 %. The high value of $RRMSE$ was associated with the large scatter verified for different strain amplitudes and the low fitting of $b_{\infty }$ throughout the test.

From the final comparison between the fitting of the constitutive model in transient conditions was verified that the non-Masing material behaviour does not allow the lower branch of the stress of the hysteresis loops to be represented with accuracy. The average value for rate saturation, $b_{\infty }$, may be used as a conservative value for the maximum and minimum stresses with deviation around 10 and 20%. The area of the hysteresis loops obtained at the end of each test is slightly larger than the experimentally verified and therefore conservative. For lower strain amplitudes and close to the elastic limit, the plasticity model presents certain difficulties in representing the shape of stress-strain hysteresis loops. Furthermore, the use of the parameter $R_{\infty ,max}$ fails to represent the effect of isotropic softening since it is twice as high as that observed experimentally. Under these conditions, the offset can be significant.

In summary, in the investigation carried out, it was possible to determine the kinematic and isotropic hardening parameters based on the Chaboche plasticity model with good accuracy for the highest strain amplitude levels. For low strain amplitude levels, close to the yield limit zone, it is recommended to use the model with caution due to the observed offsets in the comparison analysis, although these offsets make the model more conservative.