Numerical Study on the Variability of Plastic CTOD

This paper presents a numerical study on the influence of material parameters and loading variability in the plastic crack tip opening displacement (CTOD) results. For this purpose, AA7050-T6 was selected as reference material and a middle-cracked tension specimen geometry was considered. The studied input parameters were the Young’s modulus, Poisson’s ratio, isotropic and kinematic hardening parameters and the maximum and minimum applied loads. The variability of the input parameters follows a Gaussian distribution. First, screening design-of-experiments were performed to identify the most influential parameters. Two types of screening designs were considered: one-factor-at-a-time and fractional factorial designs. Three analysis criteria were adopted, based on: main effect, index of influence and analysis of variance. Afterwards, metamodels were constructed to establish relationships between the most influential parameters and the plastic crack tip opening displacement (CTOD) range, based on two types of designs: Face-Centered Central Composite Design and Box-Behnken design. Finally, the metamodels were validated, enabling the expeditious evaluation of the variability in the plastic CTOD range; in addition, the variability in the fatigue crack growth rate was also evaluated.


Introduction
Fatigue is the main failure mechanism in components submitted to cyclic loads. Design against fatigue based on damage tolerance approach assumes the presence of intrinsic defects, such as cracks. In fact, technological processes like casting, machining, welding and additive manufacturing are known to produce small defects. Accurate assessments of fatigue crack growth (FCG) rates are needed to define the time between inspections. Several models have been proposed in the literature to quantify the FCG rate as a function of loading and/or material parameters (e.g., [1][2][3]). These models are deterministic, i.e., assume that there is a well-defined relation between these parameters and FCG rate. However, the variability of the material, loading, geometry, temperature, etc., introduces uncertainty in the results. In this context, there are two main procedures to overcome the influence of sources of uncertainties [4]: (i) the use of safety factors, which is the conservative solution usually adopted to accommodate uncertainties, and (ii) performing stochastic analysis, which is often supported by probabilistic fatigue models (e.g., [5,6]), metamodeling techniques (e.g., [7][8][9][10][11][12][13]) and more recently by digital twin approaches (e.g., [14,15]). In the context of stochastic fatigue analysis, a probabilistic fracture mechanics approach to predict the fatigue life of aircraft wing attachment bulkheads was proposed by White [16]. A sensitivity analysis was performed to quantify the effect of parameters that significantly influence crack growth in metallic airframes. The studied parameters included: equivalent pre-crack size, crack growth rate, fracture toughness of the material and maximum load per flight. The variability in the total life is most sensitive to the variability in the equivalent pre-crack size, and least sensitive to the variability in the fracture toughness. In other works, FCG tests were performed on cracked specimens made of 2024-T3 aluminum alloy, illustrating the scatter in the FCG rate due to material inhomogeneity [17,18]. A predictive procedure was developed by Bahloul et al. [19], using the finite element method coupled with Monte Carlo, to evaluate the residual FCG life of a 7075-T6 aluminum alloy cracked lug component, considering the stochastic behaviour of the material parameters and the crack tip stress field. The proposed procedure showed a good ability in improving the deterministic FCG life evaluation. An uncertainty quantification methodology was proposed by Sankararaman et al. [7], to predict the crack growth as a function of the number of load cycles in a cylindrical structural component, with complex geometry and subjected to multi-axial loading with variable amplitude. The proposed methodology, which is supported by a Gaussian process metamodel trained via finite element analysis, takes into account sources of uncertainty such as loading conditions and material variability. Wentao et al. [8] proposed an approach to assess the probabilistic life of mixed mode FCG in mechanical components, by coupling finite element analysis and a Kriging-based Monte Carlo metamodel. The use of a Kriging-based metamodel for the probabilistic FCG analysis is shown to be more efficient than the Response Surface Methodology. A probabilistic fatigue life prediction model for multiple site damage in structural components was developed by Kim et al. [9], which resort to an extended finite element method coupled with a Gaussian Process metamodel for computational efficiency. Metamodels were developed by Vélez et al. [10], to relate the crack tip opening displacement (CTOD) with material properties in high-thickness offshore steel welded joints. The studied material properties (hardness, chemical composition, toughness and microstructural morphology) were shown to be significantly related with the CTOD, also being simpler and cheaper to measure and more available than the CTOD.
The aim of this paper is to carry out a numerical analysis on the influence of material parameters and loading variability in the variability of the plastic CTOD results obtained from FCG tests on middle-cracked tension (M(T)) specimens. Antunes et al. [20] developed a material law relating FCG rate with the plastic CTOD (δ p ). The FCG rate was experimentally measured in standard M(T) samples, while δ p was numerically predicted in models that simulate the experimental work regarding material properties, sample geometry and imposed loadings. It is required a certain amount of crack propagation in the numerical models for each crack length studied, in order to stabilize the cyclic plastic deformation and the crack closure level. The plastic CTOD is measured only at the end of crack propagation, and is not expected to be affected by the crack growth rate imposed artificially in the numerical model. This approach assumes that plastic deformation at the crack tip is the main mechanism behind crack tip and that the plastic CTOD is the crack driving force. Other authors proposed the use of total CTOD [21,22], but Vasco-Olmo et al. [23] proved that the plastic CTOD is the crack driving force. The da/dN-δ p model is supposed to be a material law, which can be used to estimate da/dN in other situations. In this work, the numerical models were changed regarding the material and loading parameters, and the da/dN-δ p model was used to estimate the FCG rate.
Following this introduction, the paper is organized as follows: Section 2 describes the numerical model of the FCG tests; Section 3 presents the results and discussion. Firstly, sensitivity analyses are performed to identify the most influential parameters on the plastic CTOD results. Afterwards, metamodels are constructed with the most influential parameters, to expeditiously evaluate the variability of the plastic CTOD. Finally, the proposed metamodels are validated and compared. Section 4 presents the conclusions.

Numerical Model
Numerical simulations of FCG tests were performed on a standard middle-cracked tension (M(T)) sample geometry, as shown in Figure 1. A straight crack was defined, with an initial size a o equal to 5 mm (a o /W = 0.083). Only 1/8 of the low cycle fatigue test was simulated, due to sample geometry symmetries. Plane strain state was modeled as illustrated in Figure 1b. The finite element (FE) model of the M(T) sample is composed of 6639 trilinear solid elements, with 13586 nodes. A mesh refinement was performed near the crack tip, using 8 × 8 μm 2 elements. Crack propagation was simulated by successive debonding of nodes at minimum imposed load. Each crack increment corresponds to one FE and two load cycles were applied before each crack propagation. In each cycle, the crack propagated uniformly over the thickness by releasing both current crack front nodes. 200 load cycles (i.e., 100 crack propagations) were performed, which corresponds to a total crack propagation of Δa = (100 − 1) × 8 μm = 792 μm. The first two load cycles were applied without crack increment, i.e., at a = 5 mm. The aim was the stabilization of crack opening values, since a transient behavior is observed as the residual plastic wake is formed. Additionally, the residual plastic wake may produce crack closure, which affects the crack tip fields. Again, the plastic CTOD is measured at the end of crack propagation and is not expected to be affected by the artificial crack growth rate imposed in the numerical model. The applied remote stresses σ were obtained by dividing the imposed loads by the cross-section area A (= 30 × 0.1 mm 2 ).
The numerical simulations were performed with DD3IMP in-house FE solver [24]. The average simulation time for the FCG tests is about 18 h (Intel® Core™i9-7900X 10-Core processor @ 3.3 GHz). The material elastic behavior is considered isotropic and is described by the generalized Hooke's law, where the Young's modulus (E) and the Poisson's ratio (ν) are the elastic parameters. The plastic behavior is described by the von Mises yield criterion, coupled with the Armstrong-Frederick kinematic hardening law. The von Mises yield criterion is: where σ′ is the deviatoric Cauchy stress tensor, X′ is the deviatoric backstress tensor and 0 Y is the initial yield stress. The Armstrong-Frederick kinematic hardening law is described as follows [25]: where X′  is the backstress rate, representing the translational rate of the center of the yield surface due to plastic deformation, σ is the equivalent stress and p ε  is the equivalent plastic strain rate; and the material parameters X C and Sat X represent, respectively, the saturation rate and saturation value of the kinematic hardening (exponential) law. The finite element (FE) model of the M(T) sample is composed of 6639 trilinear solid elements, with 13586 nodes. A mesh refinement was performed near the crack tip, using 8 × 8 µm 2 elements. Crack propagation was simulated by successive debonding of nodes at minimum imposed load. Each crack increment corresponds to one FE and two load cycles were applied before each crack propagation. In each cycle, the crack propagated uniformly over the thickness by releasing both current crack front nodes. 200 load cycles (i.e., 100 crack propagations) were performed, which corresponds to a total crack propagation of ∆a = (100 − 1) × 8 µm = 792 µm. The first two load cycles were applied without crack increment, i.e., at a = 5 mm. The aim was the stabilization of crack opening values, since a transient behavior is observed as the residual plastic wake is formed. Additionally, the residual plastic wake may produce crack closure, which affects the crack tip fields. Again, the plastic CTOD is measured at the end of crack propagation and is not expected to be affected by the artificial crack growth rate imposed in the numerical model. The applied remote stresses σ were obtained by dividing the imposed loads by the cross-section area A (= 30 × 0.1 mm 2 ).
The numerical simulations were performed with DD3IMP in-house FE solver [24]. The average simulation time for the FCG tests is about 18 h (Intel®Core™i9-7900X 10-Core processor @ 3.3 GHz). The material elastic behavior is considered isotropic and is described by the generalized Hooke's law, where the Young's modulus (E) and the Poisson's ratio (ν) are the elastic parameters. The plastic behavior is described by the von Mises yield criterion, coupled with the Armstrong-Frederick kinematic hardening law. The von Mises yield criterion is: where σ is the deviatoric Cauchy stress tensor, X is the deviatoric backstress tensor and Y 0 is the initial yield stress. The Armstrong-Frederick kinematic hardening law is described as follows [25]: .
X is the backstress rate, representing the translational rate of the center of the yield surface due to plastic deformation, σ is the equivalent stress and . ε p is the equivalent plastic strain rate; and the Materials 2020, 13, 1276 4 of 13 material parameters C X and X Sat represent, respectively, the saturation rate and saturation value of the kinematic hardening (exponential) law.
The material used in this study is a 7050-T6 aluminum alloy, whose elastoplastic cyclic behavior was modeled in a previous work [20]; the values of the material parameters E, ν, Y 0 , C X , X Sat that describe the elastoplastic behavior of AA7050-T6 were used as reference. In this context, a Gaussian distribution is used for describing the variability of the material parameters and the maximum and minimum imposed loads per cycle, F max and F min . The mean (µ) and standard deviation (SD) values of each parameter are detailed in Table 1, where the mean values of the material parameters correspond to the reference ones (see [20]); a coefficient of variation CV = 5% was assumed for each parameter, where CV = SD/µ. Table 1. Mean (µ) and standard deviation (SD) values of the studied input parameters; percentiles 2.5th (P2.5) and 97th (P97.5) are also shown.  Figure 2a presents the reference numerical simulation results of CTOD vs. applied stress, obtained from the mean values of the input parameters in Table 1 during the 200th load cycle. The CTOD is obtained at the first node at the left of the crack tip (8 µm from the crack tip, as schematized in Figure 2a). At minimum applied stress (point A) the crack is closed, and therefore the CTOD is equal to zero. The increase in the stress value opens the crack at point B. Afterwards, there is a linear region (point B to point C) associated with the material elastic behavior. Plastic deformation initiates after point C, and increases up to the maximum stress at point D. The dashed line shows the elastic CTOD (CTODe); the plastic CTOD (CTODp), is obtained by subtracting the elastic CTOD from the total CTOD. Figure 2b presents the evolution of plastic CTOD with applied stress during the 200th load cycle. Between points C and D, there is a progressive increase of CTODp up to the maximum stress. The range of CTODp, δ p , shown in Figure 2b, is assumed to control FCG; accordingly, the variability analysis presented in the subsequent sections will focus on δ p . The material used in this study is a 7050-T6 aluminum alloy, whose elastoplastic cyclic behavior was modeled in a previous work [20]; the values of the material parameters E, ν, 0 Y , X C , Sat X that describe the elastoplastic behavior of AA7050-T6 were used as reference. In this context, a Gaussian distribution is used for describing the variability of the material parameters and the maximum and minimum imposed loads per cycle, Fmax and Fmin. The mean (μ) and standard deviation (SD) values of each parameter are detailed in Table 1, where the mean values of the material parameters correspond to the reference ones (see [20]); a coefficient of variation CV = 5% was assumed for each parameter, where CV = SD/μ.   Table 1 during the 200th load cycle. The CTOD is obtained at the first node at the left of the crack tip (8 μm from the crack tip, as schematized in Figure 2a). At minimum applied stress (point A) the crack is closed, and therefore the CTOD is equal to zero. The increase in the stress value opens the crack at point B. Afterwards, there is a linear region (point B to point C) associated with the material elastic behavior. Plastic deformation initiates after point C, and increases up to the maximum stress at point D. The dashed line shows the elastic CTOD (CTODe); the plastic CTOD (CTODp), is obtained by subtracting the elastic CTOD from the total CTOD. Figure 2b presents the evolution of plastic CTOD with applied stress during the 200th load cycle. Between points C and D, there is a progressive increase of CTODp up to the maximum stress. The range of CTODp, δp, shown in Figure 2b, is assumed to control FCG; accordingly, the variability analysis presented in the subsequent sections will focus on δp.

Sensitivity Analysis
The sensitivity analysis enables excluding the least influential input parameters on the δ p results, to reduce computational effort in the subsequent metamodeling step. For this purpose, numerical simulations were carried out from two types of screening design of experiments (DOE): One-factor-at-a-time (OFAT) and fractional factorial design (FFD) approaches [26]. Three levels were assumed for the input parameters: the mean value (µ), the 2.5th (P2.5) and the 97.5th (P97.5) percentiles, which correspond to a 95% confidence interval. Tables 2 and 3 show the δ p results obtained from the numerical simulations under the OFAT and FFD screening approaches, respectively. Table 2. Numerical simulation results of δ p for the one-factor-at-a-time (OFAT) screening design. The influence of each input parameter on the plastic CTOD range was quantified using three analysis criteria [26]: (i) Main Effect, (ii) Index of Influence and (iii) Analysis of Variance (ANOVA). This enables assessing the influence of the analysis criterion and screening DOE on the input parameters sensitivity. The Main Effect of a given input parameter is given by: where s is the total amount of numerical simulations of the screening DOE under analysis, having input parameter values at levels P2.5 and/or P97.5 (i.e., (s = 14 for OFAT and s = 16 for FFD); δ P97.5 p and δ P2.5 p are the sum of plastic CTOD range values obtained at the 97.5 th and the 2.5 th percentiles, respectively. According to the Main Effect criterion, the input parameters are assumed influential if the absolute value of their main effect is higher than the average main effect of the seven input parameters, for each DOE.
The Index of Influence of a given parameter is given by: where δ N p is the nominal value of the plastic CTOD range, and depends on the type of DOE: In case of OFAT, it is obtained from the Simulation 1 in Table 2 (i.e., when considering the mean values of the input parameters shown in Table 1); in case of FFD, it is equal to the average of the δ p values from Table 3. According to the Index of Influence criterion, the input parameters are considered influential when the value of their Index of Influence is higher than the average Index of Influence of the seven input parameters, for each DOE.
ANOVA is used to check the significance of each input parameter, by determining the p-value for the F-test at a 95% confidence interval. ANOVA starts with the main effect of each input parameter (Equation (4)), for determining the corresponding sum of squares, SS: In the next step, the mean of squares, MS, is calculated for each input parameter: where DF is the number of degrees of freedom per input parameter, DF = l − 1, where l is the number of levels of variation considered for each input parameter (excluding the mean, µ). Since the number of levels considered is l = 2 (i.e., P2.5 and P97.5), then DF = 1. The sum of squares for error, SSE, is: where δ p i is the plastic CTOD range obtained from DOE simulation i andδ p i is the corresponding predicted plastic CTOD range value, obtained by the following linear relationship: where α j are coefficients, with j = 0, . . . , 7. The values of the coefficients were obtained by minimization of Equation (7) resorting to the Generalized Reduced Gradient non-linear optimization algorithm [27]. The mean squared error MSE is given by: where k is the number of input parameters under analysis (k = 7). The F-ratio is given by: The F-ratio enables obtaining the p-value at a 95% confidence interval. For the 95% confidence interval, the input parameters are assumed influential if the p-value is less than 0.05.
The sensitivity analysis results for OFAT and FFD are respectively indicated in Tables 4 and 5. Whatever the analysis (Main Effect, Index of Influence and ANOVA) and type of screening DOE (OFAT and FFD), the input parameters E, Y 0 and F max are consistently shown to be the most influential in the δ p results. Table 4. Main Effect, Index of Influence and ANOVA results obtained from the OFAT simulations (see Table 2). The values in bold indicate that the respective parameters are assumed influential.  Table 3). The values in bold indicate that the respective parameters are assumed influential.

Metamodeling
The variability in the δ p results was evaluated resorting to Response Surface Methodology (RSM) metamodels, constructed with the most influential input parameters, E, Y 0 and F max , to reduce the computational effort. In this context, a quadratic polynomial model is adopted to establish the following relationship between E, Y 0 , F max and δ p : where δ p RSM i is the value of plastic CTOD range predicted by the RSM model for RSM simulation i and β j are fitting RSM coefficients, with j = 0, . . . , 9. Equation (11) can be also written as a system of linear equations, as follows: where δ p i (i = 1, . . . , n) is the vector of plastic CTOD range measurements obtained from n RSM simulations, ε i is an error term and H ij is the linear system matrix, RSM simulations were performed based on two types of designs: Face-Centered Central Composite Design (FCCCD) and Box-Behnken design (BBD); the aim is to assess the influence of the choice of RSM metamodel on predicting the variability in the δ p results. Figure 3 shows the map of points (i.e., sets of input parameters) considered in the FCCCD and BBD design space, each point corresponding to one numerical simulation. Three levels of variation are considered for each input parameter: P2.5, µ and P97.5 (see Table 1); these levels are referred in Figure 3 as −1, 0, 1, respectively. The FCCCD design uses vertex points to assess the plastic CTOD range under extreme combinations between the three input parameters; as such, the resulting response surface is expected to describe extreme (and highly unlikely) behavior more accurately than that obtained using Box-Behnken (which lacks vertex points, having edge points instead). On the other hand, the Box-Behnken design is a computationally cheaper alternative to FCCCD, requiring fewer numerical simulations. The RSM numerical simulations assumed the mean values (µ) concerning the least influential parameters (ν, C X , X Sat and F min -see Table 1).
Materials 2020, 13, x FOR PEER REVIEW 8 of 13 between the three input parameters; as such, the resulting response surface is expected to describe extreme (and highly unlikely) behavior more accurately than that obtained using Box-Behnken (which lacks vertex points, having edge points instead). On the other hand, the Box-Behnken design is a computationally cheaper alternative to FCCCD, requiring fewer numerical simulations. The RSM numerical simulations assumed the mean values (μ) concerning the least influential parameters (ν, CX, XSat and Fmin-see Table 1).  Table 6 presents the least-squares solution for the RSM coefficients associated with the FCCCD and BBD designs. Whatever the design, FCCCD or BBD, the number of numerical simulations n is greater than the number of RSM coefficients β (i.e., n > j); in this case, the least-squares solution is the one that minimizes the sum of squared error at the design points,   Table 7 and Table 8 compare the δp responses obtained from the FCCCD or BBD simulations with those predicted by the corresponding RSM metamodels; the metrics R 2 and relative root mean squared error (RRMSE) were chosen to quantify the fitting performance of the design points. In general, the obtained metamodels are able to adequately describe the AA7050-T6 plastic CTOD range   Table 6 presents the least-squares solution for the RSM coefficients associated with the FCCCD and BBD designs. Whatever the design, FCCCD or BBD, the number of numerical simulations n is greater than the number of RSM coefficients β (i.e., n > j); in this case, the least-squares solution is the one that minimizes the sum of squared error at the design points, ε i = (δ p i − H ij β j ) 2 , given by [28]: Table 6. Least-squares solutions of RSM coefficients obtained from the FCCCD and BBD designs.  Tables 7 and 8 compare the δ p responses obtained from the FCCCD or BBD simulations with those predicted by the corresponding RSM metamodels; the metrics R 2 and relative root mean squared error (RRMSE) were chosen to quantify the fitting performance of the design points. In general, the obtained Materials 2020, 13, 1276 9 of 13 metamodels are able to adequately describe the AA7050-T6 plastic CTOD range measurements at the design points; also, the metamodel accuracy is best when considering the Box-Behnken design (R 2 = 0.9996, RRMSE = 0.4%) than FCCCD (R 2 = 0.9927, RRMSE = 2.0%). Table 7. Comparison between δ p results obtained from the FCCCD simulations with those predicted by the RSM metamodel; the respective R 2 value is also shown.

Metamodel Validation
The R 2 analysis shown in Tables 7 and 8 enabled a first assessment of the predictive ability of the RSM metamodels. At this stage, 60 random numerical simulations are performed to check if the response surfaces still represent a proper approximation of the plastic CTOD range within the range of variation of the main input parameters. Figure 4 presents the correspondence between numerical simulation and response surface δ p results, for the FCCCD and BBD based metamodels. In this figure, each point corresponds to one numerical simulation. The scatter of the random points around the straight line is assessed with R 2 ; the straight line indicates the ideal metamodel, in which the predicted response is equal to that obtained by numerical simulation. In general, both FCCCD and BBD metamodeling approaches provide accurate predictions of δ p , with RRMSE equal to 1.7% (FCCCD) and 0.9% (BBD), although dispersion is more evident for FCCCD (with R 2 = 0.9844-see Figure 4a) than for BBD (with R 2 = 0.9963-see Figure 4b). This corroborates the adequacy of the proposed RSM metamodels for establishing the relationship between the most influential inputs and δ p .
proposed RSM metamodels for establishing the relationship between the most influential inputs and δp.

Variability Analysis based on Metamodels
The RSM metamodels are now used to expeditiously predict the variability in the δp results. Accordingly, 100,000 random sets of input parameters E, Y0 and Fmax were generated. The values of each set of input parameters were used as input for the RSM models (i.e., see Equation (11), with the values of the RSM coefficients shown in Table 6), to predict the δp variability. Figure 5 presents the δp variability results predicted by the RSM metamodels based on FCCCD and BBD designs. Both histograms present the predicted mean δp values as well as percentiles P2.5 and P97.5, enabling the definition of a 95% confidence interval. According to the figure, the RSM metamodels predict the variability in the plastic CTOD range in a similar way, both presenting right-skewed distributions. Such skewed distributions contrast with the normal distribution, which was assumed for describing the variability of the input parameters.

Variability Analysis Based on Metamodels
The RSM metamodels are now used to expeditiously predict the variability in the δ p results. Accordingly, 100,000 random sets of input parameters E, Y 0 and F max were generated. The values of each set of input parameters were used as input for the RSM models (i.e., see Equation (11), with the values of the RSM coefficients shown in Table 6), to predict the δ p variability. Figure 5 presents the δ p variability results predicted by the RSM metamodels based on FCCCD and BBD designs. Both histograms present the predicted mean δ p values as well as percentiles P2.5 and P97.5, enabling the definition of a 95% confidence interval. According to the figure, the RSM metamodels predict the variability in the plastic CTOD range in a similar way, both presenting right-skewed distributions. Such skewed distributions contrast with the normal distribution, which was assumed for describing the variability of the input parameters.  Table 9 presents descriptive statistics of the δp histograms shown in Figure 5. The predicted mean δp value (approx. 0.340 μm) is close to the δp value obtained from the reference simulation (=0.334 μm, see Simulation 1 in Table 1), which is assumed to be the exact solution under a deterministic analysis. However, from the perspective of the current stochastic analysis, it is possible to estimate  Table 9 presents descriptive statistics of the δ p histograms shown in Figure 5. The predicted mean δ p value (approx. 0.340 µm) is close to the δ p value obtained from the reference simulation (=0.334 µm, see Simulation 1 in Table 1), which is assumed to be the exact solution under a deterministic analysis. However, from the perspective of the current stochastic analysis, it is possible to estimate not only the mean δ p value, but also the associated variance; accordingly, the predicted variability in the δ p results follows a coefficient of variation close to 15% (note that the variability in the input parameters follows σ/µ = 5%). This emphasizes the advantage and the importance of stochastic analysis, in detriment of deterministic analysis, which is key in the robust design and optimization for fatigue life of critical components subjected to time-varying loading. Table 9. Descriptive statistics of the δ p histograms shown in Figure 5.  Additionally, the following linear relationship between plastic CTOD range (δ p , in µm) and FCG rate (da/dN, in µm/cycle) from a previous work by Antunes et al. [20] allows predicting the variability in terms of da/dN: da/dN = 0.5246 × δ p (15) According to Equation (15), the predicted mean value for da/dN is close to 0.178 µm/cycle, with a coefficient of variation close to 15%. It should be highlighted that the variability analysis results may depend on the assumptions considered in the numerical model, such as the specimen geometry, number of loading cycles between crack propagations, boundary conditions, type of material, among others. This motivates further numerical studies on plastic CTOD range variability. However, linear elastic fracture parameters are commonly adopted to predict da/dN, although fatigue crack propagation is closely linked with plastic deformation occurring at the crack tip. This motivates conducting variability studies focused on the stress intensity factor range ∆K, for comparison purposes with δ p .

Conclusions
This paper concerns a numerical study on the influence of material parameters and loading variability in the plastic CTOD range of AA7050-T6 on M(T) specimens. First, the most influential input parameters were identified based on the analysis of two types of screening DOE approaches (OFAT and FFD), according to three analysis criteria (Main Effect, Index of Influence and ANOVA). Afterwards, RSM metamodels were established to predict the variability in the δ p results in an expeditious way, based on two types of designs (FCCCD and BBD). Additionally, the variability in the FCG rate, da/dN, was also predicted, based on da/dN-δ p relationships previously established in a previous work by the authors. The following remarks can be drawn:

•
The type of screening DOE and analysis does not interfere with the identification of the relevant parameters influencing the plastic CTOD range: the parameters E, Y 0 and F max are consistently shown to be the most influential; • Both FCCCD and BBD metamodeling approaches provide similar and accurate predictions of the plastic CTOD range, with RRMSE = 1.7% (FCCCD) and RRMSE = 0.9% (BBD); • The predicted variability in the plastic CTOD range results presents right-skewed distributions that follow a coefficient of variation close to 15%.