A Numerical Study of the Effect of Isotropic Hardening Parameters on Mode I Fatigue Crack Growth

The consideration of plastic crack tip opening displacement (CTOD, δp), as a crack driving force has given us the opportunity to predict fatigue crack growth (FCG) rate numerically, and, therefore, to develop parametric studies focused on the effect of loading, geometrical, and material parameters. The objective here is to study the effect of the isotropic hardening parameters of the Voce law on FCG, which are the isotropic saturation stress, YSat, and the isotropic saturation rate, CY. The increase of these hardening parameters causes δp to decrease. However, this effect is much more pronounced for YSat than CY. The variation is non-linear, and the rate of variation decreases with the increase of isotropic parameters. The increase of YSat increases the crack closure phenomenon. Finally, the influence of the isotropic parameters is more relevant for pure isotropic hardening than for mixed hardening.


Introduction
Most fatigue crack growth (FCG) studies deal with the effect of load parameters-namely ΔK, stress ratio, Kmax, and variable amplitude loading. However, the material parameters are also expected to have a major effect on FCG. In real materials, it is not possible to change material properties one-by-one; therefore, analytical or numerical approaches are the solution to study the effect of properties on FCG rate. Table 1 presents some analytical models proposed in the literature 1-10 involving material parameters-namely the elastic parameters (Young's modulus, E, and Poisson's ratio, ν), yield stress, Y0, and isotropic hardening (hardening exponent, n, cyclic strain hardening coefficient, K', and cyclic strain hardening exponent, n'). Young's modulus and Y0 are included in all of the models. The fracture toughness is included in the models of Nicholls [2] and Chand and Garg [5], while the failure strain is included in other models [3,4]. The isotropic hardening parameters are included in the models of Schwalbe [3] and Chand and Garg [5]. The effect of kinematic hardening parameters is not found in literature models. Carpinteri [10] performed finite element analyses of crack propagation by means of a strain-based criterion and observed a greater amount of crack extension in the kinematic hardening case for a given remote load, which means that isotropic hardening increases crack growth resistance. Martínez-Pañeda and Fleck [11] modeled advances in a mode I crack under small scale yielding conditions using a cohesive zone formulation. Kinematic hardening significantly raised the level of plastic dissipation. Pommier and Bompard [12] investigated the relative importance of kinematic versus isotropic hardening for modeling crack closure. The material with prevailing kinematic hardening displayed cyclic plasticity at the tip and less crack tip closure than the material with isotropic hardening. The isotropic hardening is often neglected if a material does not exhibit significant cyclic hardening and only kinematic hardening is considered [13][14][15]. In fact, studies that assume purely isotropic behavior are relatively rare [16][17][18].
The main objective here is to study the effect of isotropic hardening parameters on FCG rate. The parameters studied were isotropic saturation rate, CY, and isotropic saturation stress, Ysat. A numerical analysis based on the finite element method was followed to calculate the plastic crack tip opening displacement (CTOD), δP, which was found to be linearly related to FCG rate in numerical [19,20] and experimental [21] studies. This approach has already been used in previous works to predict the effect of a material's yield stress, Y0 [22], and Young's modulus [23]. All other physical and numerical parameters were kept constant in order to isolate the effect of isotropic hardening parameters.

Elastic-Plastic Material Model
The material studied (304L stainless steel) was assumed to have an elastic-plastic behavior. The elastic behavior was isotropic, with E = 210 GPa and ν = 0.30 as parameters of the generalized Hooke's law. The plastic behaviour was described by the von Mises yield criterion coupled with a mixed isotropic-kinematic hardening under an associated flow rule: where ′ is the deviatoric Cauchy stress tensor, ′ is the deviatoric backstress tensor described by the Armstrong-Frederick kinematic hardening law [24], and Y is the flow stress described by the Voce isotropic hardening law [25]. The non-linear (exponential) kinematic hardening model proposed by Armstrong and Frederick can be written as: where ′ is the backstress rate, which represents the translational velocity of the yield surface centre during plastic deformation, ̅ is the equivalent stress, and ̄ is the equivalent plastic strain rate; X and Sat are the kinematic hardening parameters, respectively representing the saturation rate and the saturation value of the exponential kinematic hardening X, which can be written as [26]: The Voce law is an isotropic hardening model that is often used to describe the behavior of materials that exhibit stress saturation at large strains, and can be written as follows: where YSat is the saturation stress, CY is the stress saturation rate, and ̄ is the equivalent plastic strain.
The set of material properties that best describe the cyclic behaviour of the SS304L is presented in Table 2, labeled as "Reference". The optimization procedure that was performed to identify these material properties can be found in a previous work [27,28]. Table 2 also presents the numerical changes (−25%, +25%, and +50% of the "Reference" case) in the isotropic parameters YSat and CY. The case "0.50YSat" was not considered because it leads to softening (i.e., YSat < Y0). Figure 1 presents the isotropic and mixed hardening stress-strain curves for the materials presented in Table 2, for monotonic uniaxial tension; accordingly, the isotropic hardening curves (see Figure 1a) were obtained from Equation (4) and the mixed hardening curves (see Figure 1b) were obtained from Equation (1). As expected, the level of the mixed-hardening stress-strain curves (see Figure 1b) was higher than that of the isotropic hardening curves (see Figure 1a). Accordingly, the mixed-hardening cases from Table 1 offer greater resistance to plastic deformation than the isotropic hardening cases. Moreover, the sensitivity to numerical changes in YSat is higher for isotropic hardening than for mixed-hardening (see Figure 1). Similar trends were obtained for CY.

Numerical Model
A compact tension specimen, C(T), was numerically modeled in DD3IMP using in-house code [29,30], with a width, W, equal to 50 mm, and an initial crack length, a0, of 24 mm. The symmetry conditions of the specimen's geometry allowed for the modeling of only one-quarter of the specimen, reducing the numerical overhead. Frictionless contact was modeled over a symmetry plane placed behind the growing crack front in order to simulate plasticity-induced crack closure. Relative to the specimen's thickness, t, only 0.1 mm were simulated to reduce numerical effort and to simulate a plane stress state with the proper boundary conditions. A remote cyclic load was applied at the hole of the specimen with a magnitude varying between 41.67 N and 4.167 N. The stress ratio, R, was therefore equal to 0.1 and Kmax, Kmin, and ΔK were 18.3, 1.83, and 16.5 MPa.m 0.5 , respectively. All simulations were performed with two load cycles between crack increments that occurred at minimum load. As previously mentioned, the numerical parameters were kept constant in order to isolate the effect of isotropic hardening.
An ultra-refined mesh was employed at the crack tip region with the element's dimensions being 8 × 8 μm 2 . The goal of employing this mesh was to accurately quantify strain gradients and local stress. The remaining volume of the specimen's geometry was discretized by a coarser mesh. The finite element mesh was constituted by a total of 7287 three-dimensional (3D) linear isoparametric elements and 14,918 nodes. Each crack increment-i.e., each node release-corresponded to the dimension of the elements in the ultra-refined region (8 μm). The simulations stopped when the total crack propagation, Δa, reached the value of 1272 μm, making 159 crack propagations. This allowed for the stabilization of crack tip fields and the crack closure level. Therefore, stable values of CTOD were able to be obtained. In some numerical simulations, the crack closure effect was removed by eliminating the contact of the crack flanks, which meant that the crack flanks may have overlapped. This is physically impossible; however, it can be numerically made, allowing for the study of crack propagation in a simpler state. All simulations were evaluated at the first node behind the crack tip (8 μm behind the tip). 1.25YSat 1.50YSat regime 2,3. Point 3 corresponds to the transition between purely elastic and elastic-plastic behaviour at the crack tip, and is characterized by a plastic crack tip opening displacement, CTODp, equal to 0.001 μm. This value is purely empirical. At this stage, the plastic deformation increases rapidly until it reaches the maximum value of F at point 4. The unloading of the specimens is similar to the loading process: first, the crack experiences purely elastic deformation in regime 4,5, followed by the elastic-plastic in regime 5,6. At point 6, the crack closes. Figure 2b plots a typical CTOD against the F curve in the simulations without the crack closure phenomenon. In this kind of curve, the crack is already open at the minimum value of F, making points 1 and 2 coincident. Regimes 1-3 and 3-4 are the purely elastic and elastic-plastic regimes during the specimen's loading, respectively. Regimes 4,5 and 5,6 are the same regimes, but during the specimen's unloading. The higher values of CTOD achieved in these simulations are due to the higher effective loads contributing to deformation in the absence of crack closure. This figure also shows the variation of CTODp with F. In the purely elastic regimes, CTODp remains constant and increases/decreases non-linearly in the elastic-plastic ones.  Figure 3a shows three curves of CTOD versus load for different values of YSat (cases "1.25YSat", "Ref", and "0.75YSat"-see Table 2). As expected, reducing the isotropic saturation stress has, as a consequence, higher values of CTOD, and, therefore, higher values of deformation at the crack tip. The increase of deformation is due to the plastic part, as shown in Figure 3b, which was built by extracting the elastic CTOD from the total CTOD presented in Figure 3a. In a uniaxial test, when YSat is reached, the increment in plastic deformation is performed at a constant value of stress (equal to YSat). Decreasing the value of YSat leads to a decrease in the amount of stress required for plastic strain (see Figure 1a); in other words, for a given stress value, the plastic strain will be higher for decreasing values of YSat. This can be seen by the ranges of deformations present in both Figure 3a,b, and also by the separation of the curves at lower values of the load. Similarly, for a given load value, the CTOD will be higher for decreasing values of CY.  with YSat and CY, respectively, for the purely isotropic and mixed hardening models of the SS304L. Both figures show a non-linear variation of δp with the isotropic hardening parameters, making it difficult to establish analytical models containing these hardening parameters. The increase of the hardening parameters causes δp to decrease with a rate that also decreases with the hardening parameters. The reduction of δp is much more pronounced for YSat than for CY, which can be expected when considering the physical meaning of both parameters. Also, for both figures, the mixed hardening model achieves less δp than the isotropic hardening model. This supports the results shown in Figure 1, where mixed-hardening cases (see Figure 1b) offer greater resistance to plastic deformation than isotropic hardening cases (see Figure 1a).   Figure 5a shows three CTOD versus F curves for the purely isotropic SS304L, in simulations with crack closure, in which the values of YSat were changed, corresponding to the cases "0.75YSat", "Ref", and "1.25YSat" in Table 2. The numerical changes made to YSat had a great influence on the range of CTOD values; in other words, the deformation obtained at the crack tip was deeply 304L_isotropic 304L_mixed affected. As was expected, increasing the value of YSat reduced the CTOD range and the total deformation due to the reduction of the plastic deformation, as shown in Figure 5b. Also, it can be seen here that the crack opening load, Fopen, and the crack closure load, Fclose, were affected by the numerical changes. Fopen tended to increase and Fclose to decrease when the value of YSat was higher.

Effect of YSat and CY with Crack Closure
where Fmax and Fmin are, respectively, the maximum and minimum loads. At the beginning of the numerical simulation, i.e., where Δa = 0, U * is small and increasing progressively towards a stable value, the residual plastic wake is formed. The stable values are reached at about Δa = 500 μm, which corresponds to 60 crack increments. Therefore, the 160 crack increments considered in the present study are more than enough to obtain stable values of the crack opening level. Additionally, there is more crack closure for mixed hardening than for isotropic hardening. The crack closure levels are represented in Figure 7a,b: U * is a function of YSat and CY, respectively, for both the purely isotropic hardening and mixed hardening models. As previously stated, the increase of YSat caused Fopen to increase, and, therefore, it was expected that U * would increase, as shown in Figure 7a. Also, the purely isotropic hardening reduced U * in comparison with the mixed hardening. Despite the high numerical range of variation of the isotropic parameters, the variation of U * was relatively small for both hardening conditions (approximately 17% in Figure 7a and 4% in Figure 7b). The plot of δp shown in Figure 8a,b serves as a function of YSat and CY, respectively. As noted in Figure 7a,b, the purely isotropic hardening model provided less U * than the mixed hardening model. Therefore, a more effective load range was available in the purely isotropic model, contributing to more plastic deformation, as shown in the results in Figure 8a,b. As in Figure 4a

Conclusions
A numerical approach, based on plastic CTOD, was used here to study the effect of isotropic parameters, which are the isotropic saturation stress, YSat, and the isotropic saturation rate, CY. The main conclusions are: The increase of the hardening parameters causes δp to decrease. However, this effect is much more pronounced for YSat than for CY. The variation is non-linear, and the rate of variation decreases with the increase of the isotropic parameters.
The increase of YSat causes the increase of the crack closure phenomenon. The influence of the hardening parameters is more relevant for pure isotropic hardening than for mixed hardening. Funding: This research was funded by the project no.028789, financed by the European Regional Development Fund (FEDER), through the Portugal-2020 program (PT2020), under the Regional Operational Program of the Center (CENTRO-01-0145-FEDER-028789) and by the Foundation for Science and Technology IP/MCTES through national funds (PIDDAC). The authors also acknowledge the Research Center for Mechanical Engineering, Materials, and Processes (CEMMPRE).

Conflicts of Interest:
The authors declare no conflict of interest.