An Analytical Model for the Identification of the Threshold of Stress Intensity Factor Range for Crack Growth

The value of the stress intensity factor (SIF) range threshold for fatigue crack growth (FCG) depends highly on its experimental identification. The identification and application of are not well established as its determination depends on various factors including experimental, numerical, or analytical techniques used. A new analytical model which can fit the raw FCG experimental data is proposed. The analytical model proposed is suitable to fit with high accuracy the experimental data and is capable of estimating the threshold SIF range. The comparison between the threshold SIF range identified with the model proposed and those found in the literature is also discussed. identified is found to be quite accurate and consistent when compared to the literature with a maximum deviation of 5.61%. The accuracy with which the analytical model is able to fit the raw data is also briefly discussed.


Introduction
FCG threshold (Δ th ) is one of the key parameters representing material resistance to fatigue crack growth (FCG). Newman Jr. referred to the Federal Aviation Administration (FAA) by mentioning that, traditionally, threshold is used as a limit for the damage tolerance design (DTD) [1,2]. Δ th has been used over the past 40 years in numerous FCG models available in the literature [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. However, the identification of Δ th and its application in structures residual life prediction are not quite straightforward, as it varies both qualitatively and quantitatively due to various experimental, numerical, and analytical methods and corresponding assumptions used [19]. Whilst FCG curves of physically small crack and microstructurally small crack have different shapes [20], FCG can be represented by the sigmoidal curve of log(Δ ) -log( / ) for long cracks as shown in Figure 1 [21]. This figure depicts three regions: Region I, Region II, and Region III. Region I is taken as either very slow crack growth region or near-threshold region since the SIF range of the sigmoidal curve in this region asymptotically approaches Δ th . The Paris Law [22] is normally applicable to the crack growth in Region 2. There are several models available to represent the whole sigmoidal curve covering all three regions. One of the models developed by NASA and represented by Forman and Mettu [10,21] is given in where , , , and are material constants, SIF range (Δ ) = maximum SIF ( max ) − minimum SIF ( min ), = SIF at fracture, and Δ th = SIF range at threshold.
Ideally, Δ th is the value of SIF range (Δ ) below which fatigue crack will not grow [21]. However, it has been shown [23] that cracks propagate even below the largecrack threshold measured by ASTM test procedure [1,24]. Therefore, Δ th is also defined to be a value of Δ at which crack growth rate ( / ) is below 10 −10 m/cycle [1,25]. It is sometimes experimentally determined by extrapolation to / = 0 from lower tail of the sigmoidal curve of raw data when linear-linear scale is considered [19]. Residual life of a structure can highly be influenced by a variation of Δ th . As reported by Molent and Jones [18] and mentioned by Zerbst and Vormwald [19,26], a variation in the threshold SIF range of 1 MPa m 0.5 can result in a variation of about 18% in the residual life. This also provides important insights into how relevant and reasonable it is to determine Δ th accurately in DTD approach.
There are several experimental methods available to determine Δ th . These are (i) load reduction method (LRM); (ii) max constant method; (iii) far-field cyclic compression method [19].
Load reduction method is standardised by ASTM 647 [24] or ISO 12108 [27]. The load is reduced stepwise to find Δ th in a precracked specimen at a constant . In max method, the same stepwise reduction of the load range is followed but at the same time is increased by maintaining the same maximum SIF value. The far-field compression method can be divided into three submethods: (i) Compression precracking constant amplitude (CPCA) (ii) Compression precracking load reduction (CPLR) (iii) Cyclic curve method A detailed review of all these methods is given by Zerbst et al. in [19].
The Δ th values obtained with the mentioned methods can however be quite different due to the different mechanisms involved. These mechanisms are related to the plasticity induced ahead of the crack tip as well as the conditions of the fracture surfaces. Comparatively lower threshold values have been found using the far-field cyclic compression method rather than using the load reduction method [1,2,[28][29][30]. This is due to the fact that the far-field cyclic compression method is affected by the compressive yielding at the crackstarter notch and more "steady-state" constant amplitude data in near-threshold regime is achieved with this method [29]. Crack surface roughness and grain size near the crack tip also influence the overall Δ th [1,31]. In general, greater size of grains promotes roughness induced crack closure (RICC) and oxide-induced crack closure (OICC) is enhanced simultaneously [32]. The above phenomena increase the Δ th values when measured. Consequently, in LR method, crack faces can produce rough-surface or fretting debris which contributes to the early crack closure and higher Δ th . Δ th varies with mechanically short and long cracks. Linear-elastic fracture mechanics (LEFM) is normally only applicable in long cracks under small scale yielding conditions. Newman [33] has recently referred that Δ th is not valid in gigacycle fatigue region for short cracks as there is no continuous crack propagation below ( / ) = 10 −7 mm/cycle, which is smaller than one lattice spacing per cycle [19]. In general terms, it is possible to find in the literature [34] two different Δ th levels: microstructural threshold for short crack and mechanical threshold for long crack [35]. The difference is related to the advancement of a short crack at microstructural level and stable propagation of a longer crack having a plastic zone which covers several grains. Moreover, at low SIF, the FCG rate is more sensitive to microstructure, load ratio, and environment for long cracks [20]. However, there is a minimum value independent of , which can be considered as material property and for this reason is called intrinsic threshold, also known as effective or true threshold [34]. Moreover, intrinsic threshold can be increased by the increase of stiffness and strength of the material [19,36,37]. Another important effect is related to the specimen geometry. Δ th seems to be lower in M( ) specimen than in C( ) specimen for the same Δ condition [38,39]. The justification should be related to the geometrical constraint or -stress, which is found to be lower in M( ) specimen ( -stress < 0), compared to C( ) specimen ( -stress > 0) even though -stress has different effects (e.g., PICC) which might contradict this observation. However, the lowest stress triaxiality at the crack tip associated with the M( ) specimen produces a much bigger plastic zone near the crack tip than the geometry with a high level of the constraint like the C( ) specimen [40].
Considering the fact that it is difficult to separate the extrinsic threshold from the intrinsic threshold using the crack growth data [34], the focus of this paper is to develop a model which can reliably predict the overall threshold of the material under certain loading conditions. In particular, since the model makes use of the raw data generated with a given specimen geometry under certain loading conditions, the analysis of the raw data includes both the load ratio and the -stress effects. The value identified with the model can be an intrinsic or an extrinsic value depending on the testing conditions at which the data has been acquired.
As discussed above, the Δ th value usually decreases with the increase of [41]. Two types of -dependency have been reported in the literature [34]. In some cases, Δ th decreases up to a critical value of and then it becomes constant beyond that value [19]. In other cases, Δ th continues to decrease beyond the critical value of [42]. Klesnil and Lukáš [3] used the following equation to identify Δ th in a steel alloy: where is the stress ratio, Δ th0 is the fatigue threshold value at = 0, and is the material constant. However, other approaches [43] have been adopted like the one reported by Kwofie in which an equivalent stress approach based on ratio is used to identify fatigue threshold value. In general, it has been recognised that crack closure is found to be the controlling factor in this case [18,44]. For this reason, a different parameter has been introduced Δ thr , which is a FCG threshold value that depends on and the crack length value. In the literature, the scatter in fatigue life was explained by the variation of Δ thr values [18]. Further methods to experimentally identify the threshold condition have been recently developed using plain fretting crack arrest analysis. The dispersion between long crack Δ th fretting estimations and conventional fatigue data was found to be less than 10% [45].
Due to the high variability of the Δ th values, the determination of the FCG threshold cannot be certain [19]. Although several models have been proposed to experimentally identify the threshold values, all of them suffer from issues related to the plasticity induced closure effects. For this reason, threshold values reported in the literature for the same material can vary in a wide range due to the different procedures that were followed. The aim of this paper is to present a new procedure to identify FCG threshold value for long crack which can overcome the problems related to the experimental procedures reported in the literature. The analytical model proposed here makes use of FCG data obtained from -increasing tests, which are used to derive the FCG properties of the material under the long crack condition, allowing us to identify under the same testing conditions the three regions of the entire sigmoidal curve, from the threshold condition up to the final value of the crack length.

Test Results for Model Development
Propagation models built on results obtained from a limited number of tests not only have a validity range closely linked to the particular experimentation carried out, but also are not suitable to fit all crack growth data with the same accuracy for the whole field of number of cycles for each test [46]. In order to overcome these drawbacks, several FCG datasets obtained with different materials, loading conditions, and types of specimens have been collected from the literature. These datasets have been used to verify the suitability of the model in fitting the experimental raw data as well as identify the Δ th values of the materials at the corresponding values. A short description of the datasets collected from the literature is as follows. [47]. Ghonem and Dore [47] carried out tests at room temperature using M( ) specimens made  of aluminium alloy 7075-T6 having a thickness of 3.175 mm.

Ghonem and Dore
The crack direction was perpendicular to the rolling direction and the loading conditions are reported in Table 1. Sixty specimens were tested under each loading condition. [48]. The experimental activity reported by Virkler et al. [48] was aimed at determining which crack growth rate calculation method yields the least amount of error when the crack growth rate curve is integrated back to obtain the original " " versus " " curve data.  [49]. The experimental work of Wu and Ni [49] was carried out on compact tension C( ) specimens made of aluminium alloy 2024-T351, having thickness = 12 mm and width = 50 mm. Tests were carried out with variable and constant amplitude loading. The two samples marked by the authors as CA1 and CA2 and composed of 30 and 10 specimens, respectively, were tested at constant amplitude loadings reported in Table 2.

Model Implementation
The analysis of experimental data obtained from FCG test is quite complex due to the scatter nature in the raw data which is amplified by the derivation needed to compute the FCG rate. Several useful formulae to fit the experimental data with the aim of a better, smoother curve have been proposed and reported in the literature. Among those, the use of a polynomial function to fit the raw data gives the possibility of obtaining a single numerical expression of the crack growth rate valid in the entire data range [46]. The choice of the most appropriate function can be made considering that the crack growth is exponential by nature. In mathematical terms, an exponential correlation can be represented introducing logarithmic functions for the crack length [50][51][52]. This linear correlation (log( ) versus ) can be represented on a semilogarithmic plane as a straight line. There are models proposed in the literature which are developed adopting an exponential structure [52]. However, the trend identified using the experimental FCG data changes as the crack length approaches the failure condition. This consideration is supported by the presence of three different regions in the sigmoidal curve with each of them following a different trend. On the basis of the aforementioned observations, the most suitable formula to fit the whole FCG experimental data points can be deduced by summating the individual effects of the different crack growth regions [53]. Therefore, the following model, on the basis of a trial-and-error method, could be established: where , , and are three parameters to be determined by the least-square method. The procedure to derive the values corresponding to ℎ and th is described in later parts of this paper. The proposed model makes use of nondimensional fatigue crack life, which makes it more general. Moreover, the nondimensional fatigue crack life allows decoupling the identification of the equation parameters, which are meant to be a material property, from the actual total life for the particular test. The nondimensional fatigue life is defined as follows: The parameter th , which is identified through best-fit curve together with three parameters ( , , and ) reported above, is related to the nucleation phase and hence to the threshold value. is the final value of the experimental crack life, which is the number of cycles counted from the initial crack length up to the final failure of the specimen, whilst is the generic value of the fatigue crack life.
Useful formulae can be derived for other parameters in (3) by considering some specific data points of the crack growth curve. At = , which corresponds to the last experimental data point of the test, the crack length is equal to the value of the crack length in the corresponding last front just before the failure condition of the specimen. This gives Similarly, considering the value of (3) at = 0, which corresponds to the first experimental data point, the crack length is equal to the value of the crack length th corresponding to the starting point of the test. This gives As already stated, the parameter th is related to the threshold condition and represents the number of cycles needed by the crack to reach the crack length corresponding to the threshold condition. From (3), the crack length in correspondence to the threshold condition is equal to the value of the th parameter in correspondence of = − th . Equation (3) is a continuous differentiable function in the range th < < . It is therefore possible to derive the analytical expression of the crack growth rate ( / ) as a function of . The function (see (7)) can be used to represent the continuous propagation process from threshold region up to the final fast crack growth region ) ⋅ (( + th )/( + th )) /( −(( + th )/( + th )) ) .
The parameters in the crack growth rate function are identified by means of the linear regression using the FCG raw data. The analytical expression of the crack growth rate is equal to zero at = − th according to the assumption that the crack length at this value corresponds to the threshold condition. The procedure of applying the formulae of the analytical model to derive the threshold SIF range is as follows: (i) The experimental raw data of crack length versus number of cycles are fitted using (3). The method adopted for the fitting is the linear regression to identify the four parameters of 0 , , , and and minimize the error. In an earlier paper [53], the model here presented was adopted to assess the accuracy in fitting the raw data produced during FCG tests. In the same paper [53], the normal distribution of the residuals as well as the distribution of the equation parameters has been included. In the present paper, the discussion is focused on the identification of the threshold SIF range through the use of the analytical model proposed by the authors.   (7). The crack length values derived from (3) are used to deduce the SIF range values by means of the expressions in accordance with the international standards. This means that the method requires the knowledge of the closed form of the SIF for the tested specimen. In this paper, the expressions reported by the ASTM E647 [24] have been used to compute the SIF values.
(v) The value of the crack length corresponding to = − th is used to derive the value of the threshold SIF range which corresponds to a crack growth rate equal to zero.
The procedure described above has been implemented in a Matlab code to identify the FCG curves and the Δ th values using the datasets produced by Ghonem and Dore, Virkler, and Wu and Ni. The detailed discussion about the capability of the model to properly fit the datasets used in this paper is given in an earlier paper [53]. Figure 2 shows some examples of fitting results with the experimental points related to Set I produced by Ghonem and Dore.
Moreover, normality of the residuals obtained from each curve has been verified by the 2 normality tests and the corresponding residuals frequency histograms have also been evaluated. In Figure 3, the means of residuals for Ghonem and Dore Set I and Set III are shown as an example to highlight the notion that the mean value is equal to zero [53].
A further version of the Matlab code, which was already implemented for the fitting curves, was developed further in order to identify the values of the FCG rate as well as the SIF range values. In particular, the SIF values in correspondence to = − th for each curve of all datasets have been computed in order to estimate the threshold values and compare these with those reported in the literature.

Results and Analyses
The interpolation of the raw experimental data represents the first step of the analysis. The suitability of the equation for fitting the data has been summarised in the above section. In particular, raw data fitting has an average value of 2 equal to 0.9998 for all datasets [53]. The values of the parameters identified for each dataset are shown in Table 3 The values reported in Table 3 have been computed as an average of the values identified over the total number of tests for each dataset.
The crack length as a function of the number of cycles derived in the range [− th ; ] is shown for each dataset in Figure 4. In particular, the curve fitting related to the three datasets produced by Ghonem and Dore is shown in the top row of Figure 4 (Set I, Set II, and Set III), and the curves related to the dataset produced by Virkler and the curves related to the datasets produced by Wu and Ni are shown in the bottom row of Figure 4. In each dataset, all the fitting curves tend to the same asymptotic value as the number of cycles approaches − th . However, the values are different between the various datasets. In order to make the comparison between the experimental data points and the curves identified with the analytical model possible, the logarithmic plot of the horizontal axis, which is the number of cycles, has been used. The logarithmic scale is only used for the sake of clarity of the graph whilst the equations adopted are not affected by this choice.
In order to draw the FCG curve for the entire range, it is necessary to derive the crack growth rate together with the SIF range for the corresponding values. The curves shown in Figure 5 correspond to all experimental data of the datasets considered in this paper. These graphs show clearly that the gradient approaching = − th is equal to zero, which reflects the asymptotic behaviour in the -curves. As a consequence, the FCG rate, as expected, approaches zero.
This observation can be used to extrapolate the crack growth curve from the lower part. The value of the threshold SIF range is found where = − th . By means of (3), the value of the crack length at = − th can be derived. In Figure 6, the values of the threshold SIF predicted by the model for each curve of the five datasets analysed in this paper are shown together with the corresponding values gathered from the literature [44].

Discussions
The curves of crack length versus number of cycles found from the model correlate well with the raw data points. To relate the parameter values obtained by means of the best fit to the testing conditions, the contour plots shown in Figure 7 have been created. The contour plots can be used to identify the range in which the optimal values of the model parameters should be identified by means of the best fit of the experimental raw data in terms of crack length versus Ghonem and Dore Set I      number of cycles. Nevertheless, the possibility of correlating the model parameters with the loading conditions by means of analytical expression can be adopted. This would allow us to either reduce the number of parameters that can be identified through the best fit or identify the FCG curve for a given material under certain loading conditions, which are usually defined in terms of maximum load and load ratio. For this reason, the correlations between the 6 parameters of the proposed model and the maximum load max and the load ratio have been investigated using the values obtained from the fitting analysis carried out with the datasets found in the literature and already discussed in this paper. In (Figures 7  and 8), max is defined as the maximum value of the applied load in each cycle. The two parameters and and th depend on the load ratio for low values of the maximum applied load only. At higher values of the maximum applied load, the range of the above parameters is the same whatever the load ratio is, which means they are independent of the load ratio. Moreover, their values decrease as the values of and max increase. A different trend can be observed from referring to the model parameters of th , p, and ℎ. The combined effect of the load ratio and the maximum applied load is always present for the investigated range values of the loading conditions. In addition, the values of the two parameters and ℎ decrease as the values of and max increase, whilst for th the trend is opposite with regard to max . In the graphs shown in Figure 8, the correlations identified between the model parameters and the loading conditions are shown. In some cases, the combination of more than one parameter has been considered.
In particular, for ℎ, th , and , three different expressions have been introduced: Expressions reported as (8) and (10) have been shown to be related to the maximum applied load, as shown in Figure 8. Equation (9) is related to the load ratio and it is shown in Figure 8. In terms of the quality of fit, it can be assessed with the value of the coefficient of correlation ( 2 ). Good fitting is achieved if the coefficient of correlation is between 0.8 and 1 for all the correlations identified [54]. th in function of max , as shown in Figure 8, has lower coefficient of determination of 0.8286 in the current study. This value is however still within the acceptable range and the correlation between the considered parameters values is acceptable.
The extrapolation of the curves as shown in Figure 5  reaches = − th asymptotically, which can be used to generate Δ th values based on the theory of identifying threshold [19,21]. The range between the maximum and the minimum values of the predicted Δ th for each dataset is quite small so that the average value of the threshold is used. The reference values for materials considered in this paper have been taken from the literature in order to verify the model. In the literature, a range of threshold values are 8 Advances in Materials Science and Engineering  Figure 6: Threshold SIF range for the five datasets.
also mentioned for a particular material for a constant [18,41,55]. Figure 10 shows the comparison between the threshold values obtained using the model and threshold values found in the literature for the same kind of material. The values considered as literature values are based on ESDU documents [41] for aluminium alloy 7075-T6 (at = 0.4, 0.5, and 0.6) and aluminium alloy 2024-T3 (at = 0.2). For 2024-T351 (at = 0.2), the literature value is taken based on Δ th value provided by [44] and normalized threshold SIF range against curve provided by [43]. The average threshold value of the model is used for the comparison. In Ghonem and Dore (Sets I-III), Virkler, and Wu and Ni datasets, the threshold values of the corresponding materials from the literature are 1.5 MPa√m, 1.7 MPa√m, 1.9 MPa√m, 2.85 MPa√m, and 3.6 MPa√m. Comparing these values with those predicted by the proposed model, the percentages of error are approximately 5.3%, 4.1%, 0%, 2.1%, and 0.3%, respectively.
As shown earlier in Figure 6, in all cases, there is a band of threshold values predicted for the same material with the literatures showing the same trends. The threshold versus graph provided by [41] shows a small range of threshold values with upper and lower limit indicated for 2024-T3 aluminium alloy. Different investigations by Newman Jr. and Ruschau [1,2,[28][29][30] also found a range of thresholds based on different experimental methods used, for example, load reduction (LR) methods, max constant methods, and farfield cyclic compression methods.
Recently, Molent and Jones [18] used a FCG threshold parameter Δ thr to explain the scatter in the fatigue life prediction. Δ thr is used to explain the dependency of the threshold values on the material properties, ratio, crack length, and loading method used for the testing phase. From the literature [45,56], it was found that Δ thr = 0 produces conservative crack growth predictions for different aluminium alloys. The fact that different values of the threshold SIF were reported for the same material in the literature can be treated as supporting evidence to the concept of cyclic stress intensity threshold. For the material with which the samples used by Virkler were made, a range of Δ thr (2.9 ÷ 4.2 MPa√m) is reported in the literature, whilst, in the material used by Wu       Load ratio R Figure 9: Contour plot of Δ th in function of the maximum applied load and the load ratio.
= 0.33 for 2024-T351 aluminium alloy has been shown. All these values are in good agreement with those estimated by the proposed model.
In order to investigate the variability of the threshold SIF range, the response surface reporting the correlation of this latter parameter as functions of the maximum applied load and the load ratio is shown in Figure 9. The values shown in the graph related to the threshold SIF are those identified through the fitting analysis of the raw data. Although all the values identified in the analysis related to different aluminium alloys are shown on the same plot, the validity of the contour plot is limited to the and max ranges investigated in this paper and the alloys considered in the tests are referring to the datasets used for the analysis. However, the contour can represent a valid tool for a preliminary identification of the threshold SIF range.
Reporting all the results on the same graph in Figure 10, it is possible to identify the common trend useful to compare the results with the literature values. Firstly, the threshold line found with declined linear pattern or shape in relation to for the 7075-T6 aluminium alloy is qualitatively and quantitatively consistent with the line found in the literature [41]. The percentage of error ranges between 0.24% and 5.61%, which is quite low considering the scattering nature of the fatigue test data. The threshold SIF range values should converge at the higher value of but it was found that the scatter was getting bigger and reached 5.61% at = 0.6. This difference or scatter could be explained by different experimental methods used and the corresponding crack closure effects as referred by [1,2,[28][29][30]. All the predicted threshold values of aluminium alloys (7075-T6 and 2024-T6) underestimate the values from the literature (see, e.g., [41]) except for the threshold value of 2024-T351 aluminium alloy which overestimates the literature value from [43,44]  derived using the analytical model here proposed is equal to 2.3651 MPa√m with an error less than 3%.
In general, the yield strength of 7075-T6, 2024-T3, and 2024-T351 aluminium alloys is 510 MPa, 350 MPa, and 330 MPa, respectively. The threshold values found from the model for the three alloys at = 0.2 are 3.6 MPa√m, 2.8 MPa√m, and 2.4 MPa√m, respectively. This indicates a correlation between the yield strength and threshold of these aluminium alloys: higher strength aluminium alloy possesses relatively higher threshold and vice versa. This supports the statement of the previous investigations which found the same type of correlation between the strength and threshold of materials [19,36,37].
It should be added that the values of th in (6) were found to influence the threshold value since it is the value of the crack length corresponding to the fatigue life value equal to = − th . Therefore, th values can be further considered to correlate with the material properties and the load ratio as well as the geometry of the cracked component. It is possible to find in the literature that the threshold value increases with the strength of the material [19,36,37]. Further research should be carried out to properly address this correlation, which could help in both identifying the parameter values and giving them a more physical meaning.

Conclusions
An analytical model for the interpolation of crack propagation data has been developed. The threshold SIF range has been derived for different materials and different specimen geometries. The new model has been shown to fit, with the needed accuracy, to a wide range of experimental data produced with different specimen geometries, different materials, and different loading conditions. Moreover, it has been highlighted that it is possible to identify, by means of the above model, the value of the threshold SIF range with an error, as compared to the values reported in the literature, of less than 6%. The relation between Δ th and ratio predicted by the model agrees well with the literature results. The proposed model can therefore be valuable in identifying the threshold of stress intensity factor range for fatigue crack growth.

Nomenclature
: Crack length (mm) : Final value of the experimental crack length (mm) : Specimen thickness (mm) , , , : Forman's constants / : Crack growth rate (mm/cycle) ℎ, , , th , , , th : Parameters of the model proposed in the present paper : F r a c t u r e t o u g h n e s s ( M P a √mm) Coefficient of determination : Specimen width (mm) : Material parameter for Klesnil and Lukáš's model Δ : Range of the stress intensity factor (MPa√mm) Δ th : Value of the stress intensity factor range threshold (MPa√mm) Δ th0 : Value of the stress intensity factor range threshold for = 0 (MPa√mm).