**3. Experiment**

Tubular unnotched thin-walled specimens made of CuZn37 brass were subjected to fully reversed constant-amplitude fatigue loading under various loading paths. Because the experimental test details can be found in [64], only the most important information is reported. The strain-controlled fatigue tests were conducted according to the ASTM E2207-02 standard. The failure definition was a 10% drop in the axial force or torque. The monotonic and cyclic mechanical properties of the CuZn37 brass are listed in Table 2. The regression fatigue curves (13) and (14) are shown in Figure 1. The Monte Carlo sampling technique was applied to generate the distribution of the fatigue curves ((*i*)-indexed curves in Figure 1) using a multivariate normal distribution of regression coefficients [36,65]. A sample size of 1000 was applied, estimated from the stability analysis of error indexes (Section 4.1). The blue vertical dashed lines in Figure 1 indicate the overlapping experimental fatigue life regimes of the axial and shear fatigue curves. The experimental setup included eight different loading paths with a fixed strain ratio of amplitudes *<sup>γ</sup>xy*,*<sup>a</sup>*/*<sup>ε</sup>xx*,*a*. Each loading path was tested using 10 specimens at different strain loading amplitudes. The shapes of the loading paths with applied strain ratios are presented in Figure 2a. The material responses in the form of stress paths are displayed in Figure 2b. The stress amplitudes implemented in the fatigue life prediction were determined for the half-life. The received experimental fatigue lives 2*Nexp* were within the range of [383, 191000] reversals. The experimental data file with the registered signals (strain and stress components) was uploaded to the Mendeley data repository (https://data.mendeley.com/datasets/7fbkf6y7gv/1, accessed on 25 October 2021).

**Table 2.** Mechanical properties of CuZn37 brass.


**Figure 1.** Shear and axial strain-life fatigue curves for CuZn37 with error indexes *MPE*, *SD*, *T*(0.95) (described in Section 4.1).

**Figure 2.** *Cont*.

**Figure 2.** Registered experimental strain paths in (**a**) axial-shear strain space, *εxx* − *γxy* and (**b**) axialshear stress space, *σxx* − *<sup>τ</sup>xy*.

### **4. Results and Discussion**

The normal and shear stress/strain components on the critical plane implemented in the parametric (Section 2) and GP-based fatigue life prediction (Appendix A) models were calculated based on the registered strain and stress tensor components. Searching for the plane orientation with the maximum shear strain (the critical one), numerical simulations were conducted with the application of Euler angles, *ϕ*, *θ*—describing the plane orientation [66]. In the simulation, the step between subsequent values of the Euler angles was equal to 1◦. For the non-proportional loading paths, instability in the calculated fatigue lives was detected as a result of the existence of different planes with the same maximum value of shear strain but with different normal stress and strain components. To overcome this problem, the critical plane was searched within the boundary 1 − *δ*, <sup>1</sup> max *ϕ*,*θ γns*, where *δ* = 0.0001 was applied. Within this boundary, the maximum values of normal stress, normal strain, and shear stress were determined and applied to the fatigue life calculations (the values could be found in the file uploaded to https: //data.mendeley.com/datasets/7fbkf6y7gv/1,accessedon25October2021).

### *4.1. Error Indexes*

The calculated fatigue lives were compared with the experimental fatigue lives to estimate the performance of each applied fatigue model. A few error indexes were implemented in which the first group belongs to purely statistical parameters of fitting, and the second group belongs to parameters to estimate the applicability of models to life prediction. The first group is based on a percentage error, defined as

$$PE = \frac{\log \text{N}\_{\text{exp}} - \log \text{N}\_{\text{cal}}}{\log \text{N}\_{\text{exp}}} \times 100,\tag{20}$$

where *Nexp* and *Ncal* are the experimental and calculated number of cycles to failure, respectively. The statistics of *PE* provide information on the efficiency of the model in fatigue life prediction. The mean values, *MPE* and standard deviation, *SD* of *PE* were evaluated. The positive values of *PE* concern cases where the calculated life was shorter than the experimental one (conservative estimate), whereas non-conservative estimates were indicated by negative *PE* values. An *MPE* value equal to zero indicates a perfect prediction-unreal case owing to fatigue life scatter.

The mean and standard deviation of *PE* estimate the fitting properties of the applied model. However, to build (based on these parameters) the ranking of models with respect to their acceptable uncertainty in life prediction can be questionable. First, the two parameters are combined into a decisive one. Second, the values of these parameters are not intuitive if the predicted fatigue lives are acceptable. The experimental fatigue life can vary significantly even under a single loading condition [67,68]. The fatigue scatter factor is defined as: 

$$T = \begin{cases} \begin{array}{c} \frac{N\_{exp}}{N\_{cal}} \quad for \quad N\_{exp} \ge N\_{cal} \\\ \frac{N\_{cal}}{N\_{exp}} \quad for \quad N\_{exp} < N\_{cal} \end{array} \tag{21}$$

Fatigue scatter factor less than three is commonly accepted [69,70], and a value less than or equal to two exhibits a very good estimation [71–73]. Thus, the statistical distribution of the *T* factor was used to define [74] a more intuitive error index as a 0.95-quantile of the fatigue life scatter factor. The 0.95-quantile (*T*(0.95)) was calculated based on the shape-preserving piecewise cubic interpolation of the empirical cumulative distribution of *T*. The value of *T*(0.95) determines the minimum fatigue-life scatter band required to include 95% of all experimental fatigue data.

### *4.2. The Parametric Fatigue Life Prediction Models*

A comparison of the experimental 2*Nexp* and calculated 2*Ncal* fatigue lives for the four analyzed parametric models is presented in Figure 3. Each panel in Figure 3 includes the solid line of perfect life consistency enclosed by dashed lines of the scatter band *T* = 2.0. The error bars indicate 95% prediction intervals computed using the Monte Carlo sampling technique on strain-life fatigue curves. Additionally, the results for different loading paths are marked, and the corresponding values of the error indexes are included in the legend of the figure panels. The data used for the calibration of the model, i.e., uniaxial tension– compression and pure torsion were labeled as 'Train data' and the remaining data as 'Test data'. For the labeled data, the *MPE*, *SD*, and *T*(0.95) error indexes were estimated and presented in each panel.

The error indexes are identical for the analyzed models of Fatemi–Socie, Brown–Miller, and Glinka et al., with the application of the life-dependent material parameters that equalize the performance of the parametric models under data used for calibration. The obtained error indexes *MPE* = −0.1%, *SD* = 2.5%, and *T*(0.95) = 1.5 for the training data and the models characterize the experimental life scatter for the CuZn37 brass under the fatigue test conditions. These values were treated in further analyses as the reference consistency of fatigue lives. Yu et al. applied constant weighting factors of shear and normal strain energy parameters, and the results for the analyzed CuZn37 brass were ineffective even under uniaxial and pure torsion loadings with *MPE* = −23.7% and *MPE* = −9.8%, respectively (Figure 3d).

The fatigue lives for the multiaxial proportional loading path were calculated using models with life-dependent material parameters within the acceptable value of error indexes (80% of the data is included within the scatter *T* = 2.0, and 100% of the data is included within the scatter *T* = 3.0). For the non-proportional loading paths, the best life prediction performance was exhibited using the Fatemi–Socie model (Figure 3a). However, the fatigue life under loading path NPR4 was unsuccessfully estimated with *MPE* = 7.7%. The performances of the models by Brown–Miller and Glinka et al. were the worst, with *MPE* values exceeding 14% and 27% for the NPR4 path, respectively. The Brown–Miller model with shear strain *γns* linearly combined with the normal strain *εn* resulted in a very conservative life estimation under non-proportional loading paths (Figure 3b). This means that the non-proportional loading paths led to a higher value of the normal strain on the critical plane of the maximum shear than expected using the Brown–Miller model. The same effect, but magnified by the shear stress on the critical plane, was observed in the model by Glinka et al.

**Figure 3.** Comparison of experimental and calculated fatigue lives for the models of (**a**) Fatemi–Socie, (**b**) Brown–Miller, (**c**) Glinka et al., and (**d**) Yu et al.

The tendency of higher life underestimation (*Ncal* < *Nexp*) for longer experimental fatigue lives was observed for all models with the life-dependent material parameters under non-proportional loading paths. It is concluded that the CuZn37 brass experienced additional hardening under non-proportional loading paths that resulted in higher stress values than expected by the models, which were calibrated by uniaxial and pure torsion loading paths. Yu et al. overestimated (*Ncal* > *Nexp*) the predicted fatigue lives for all loading paths.

### *4.3. The Gaussian Process-Based Fatigue Model*

### 4.3.1. Physics-Based Input Parameters (Predictors)

The GP model is assumed to substitute the semi-empirical fatigue criteria; thus, the scalar output of the logarithm of the fatigue life must be invariant under the rotation of coordinate systems. To meet this requirement, the covariance function for the GP must not directly implement the input parameters as stress and/or strain tensor components. The rotation of the coordinate system transforms the stress/strain tensor components, and thus, the model trained on the dataset valid for the original coordinate system would not apply to the rotated coordinate system.

A state-of-the-art review of fatigue life prediction models (Section 2) indicates a few commonly applied crucial physics-based quantities. These are derived from the concept of the critical plane built on the experimental observation that fatigue cracks in metals are initiated in the plane of maximum shear stress or strain [49,75], *τns*, *γns*. Furthermore, these microcracks could not be developed if the maximum normal stress or strain *σn*, *εn* on this plane were below the critical value. These four parameters are considered primary in most multiaxial fatigue life prediction models (Table 1). Owing to the zero mean stresses in the experimental data, the *<sup>σ</sup>n*,*<sup>m</sup>* predictor is omitted for components of the input data vector *x*, as follows:

$$\mathbf{x} = [\gamma\_{ns}, \ \varepsilon\_n, \ \tau\_{ns}, \sigma\_n]. \tag{22}$$

The maximum shear strain plane was determined over the entire loading history and fixed in the fatigue life calculation process. Thus, the selected predictors (22) are invariant under the rotation of the coordinate system. Regarding the identification of multiple planes with equal maximum values of shear strain, the plane with the highest normal stress was selected as the critical plane. The relevance of each predictor was analyzed during the training process, as presented in the next section. The diagram of data flow for the fatigue life prediction based on the proposed GP fatigue model is presented in Appendix A, Figure A1.

### 4.3.2. Training Process

The polycrystalline structure with preferred slip systems makes metallic materials sensitive to the type of fatigue loading paths. For example, the non-proportional loading characterized by the rotation of principal stresses could activate a larger number of slip systems and their intensive interactions compared to the proportional loading path [76]. These phenomena could influence the fatigue life of materials; however, this effect depends on many factors [77,78]. If the fatigue GP-based model is expected to be implemented under non-proportional loading for materials sensitive to its effect, the training data should include such a case. The analyzed experimental data included various non-proportional loading paths (Section 3, Figure 2), and the most common one (easy to replicate) with a 90◦ phase shift (Path NPR in Figure 1) was selected for inclusion in the training process. To analyze the effect of the non-proportional loading path on the fatigue life, two training datasets were implemented. The first training dataset includes only loading paths commonly applied to identify the parametric fatigue life prediction models (Section 2), that is, uniaxial push-pull (10 specimens) and pure torsion (10 specimens) loadings. The second training dataset additionally included the NPR path (10 specimens).

Initially, the training process included all four predictors (22), five specified covariance functions (Appendix B), and the first training dataset. The length scales *li* for each covariance function and each predictor were used to calculate the relevance factor (*RF*), defined as:

$$RF = \left(\frac{l}{std(\mathfrak{x})}\right)^{-1} \text{ .} \tag{23}$$

where *l* is the length scale, and *std*(*x*) is the standard deviation of the analyzed predictor observed in the training data. The larger the length-scale parameter, the lower the covariance, and thus the lower the influence of its predictor. However, it also depends on the analyzed physic-based quantity; for example, the strains in the fatigue regime are several orders of magnitude below the stresses. This can be considered by normalizing the length scale using the standard deviation of the input of the analyzed predictor, *std*(*x*). The reverse of *l*/*std*(*x*) is defined as the *RF*. An *RF* value approaching zero indicates that the analyzed predictor can be ignored.

The RFs obtained for the first training dataset are shown in Figure 4. The vector of the standard deviations for the first training dataset is std(*x*) = [0.0026 (−), 0.0009 (−), 6.6 (MPa), 62 (MPa)]. Two conclusions can be drawn. First, the exponential covariance function (EX) exhibits the lowest values of the relevance factors compared to other kernels. This means that the predicted values for the test data applying the EX kernel will be characterized by a larger uncertainty than for the remaining kernels. Second, the predictor *<sup>ε</sup>n*—of the normal strain on the critical plane can be neglected for M5/2, RQ, and SE kernels (near zero relevance factor). The mean percentage error (MPE), standard deviation of (PE), and *T*(0.95) factor for all the covariance functions are in the ranges *MPE* = [−0.07, −0.06] %, *SD* = [2.0, 2.5] %, *T*(0.95) = [1.53, 1.54]. Based on the obtained results, further analysis was reduced to four covariance functions (M3/2, M5/2, RQ, and SE) and three predictors, *x* = [*<sup>γ</sup>ns*, *τns*, *<sup>σ</sup>n*]. The error-fitting indicators; *MPE*, *SD*, and *T*(0.95) were in the same range as the initial selection. The estimated values of the hyperparameters for the analyzed kernels for the first training dataset are presented in Table 3.

**Figure 4.** RFs obtained by application of different kernels for the first training dataset (uniaxial and torsion).


**Table 3.** Hyperparameters of the covariance functions for the first training dataset.

The relevance factors for the second training dataset, including the non-proportional loading path NPR, are presented in Figure 5. Based on the previous analysis, estimation was conducted on four covariance functions and three predictors. Here, the vector of standard deviations of the predictors was found as std(*x*) = [0.0027(−), 22.4(MPa) <sup>124</sup>(MPa)]. For the second training dataset, all four kernels indicate the irrelevance of the shear stress *τns* on the critical plane. Thus, in the final selection, the shear stress as a predictor was abandoned. The error fitting indicators for both selections are in the ranges of *MPE* = [−0.07, −0.06] %, *SD* = [2.15, 2.50] %, *T*(0.95) = [1.44, 1.67] %.

**Figure 5.** RFs obtained by application of different kernels for the second training dataset (uniaxial, torsion, and NPR path).

Inclusion of the non-proportional loading path NPR into the training data increased the RFs for all kernels by approximately four times. The estimated values of the hyperparameters for the analyzed kernels for the second training dataset are presented in Table 4.

**Table 4.** Hyperparameters of the covariance functions for the second training dataset.


### 4.3.3. Test Process

In contrast to the parametric models in which estimation of output uncertainty requires additional methodologies, for example, Monte Carlo sampling [36], the inherent property of the GP model is the variance estimation of the outputs. This property was utilized in the fatigue life prediction shown in Figure 6 by additional vertical bars providing 95% prediction intervals. Figure 6 presents a comparison of the experimental and calculated fatigue lives obtained by the implementation of M3/2, M5/2, RQ, and SE kernels trained on the first dataset (only tension-compression and pure torsion paths).

The results demonstrated that within the test data (6 loading paths and 60 specimens), only fatigue lives predicted for the multiaxial proportional loading can be accepted with *MPE* = [3.7, 4.2]% and all data included within the fatigue scatter band *T* = 3.0. The fatigue lives under the non-proportional loading paths with no exception are predicted with high underestimation with *MPE* = [22.8, 37.6]% and *T*(0.95) = [13.2, 88.4]. Additionally, the 95% prediction intervals are approximately two and five times larger than those estimated for the training data with the implementation of the M3/2 and M5/2 kernels. A summary of the error indexes received for the analyzed kernels is presented in Figure 7. The lowest error index was obtained for the M3/2 kernel, but with high (95%) prediction intervals and unacceptable values of *T*(0.95) > 3.

**Figure 6.** Comparison of experimental and calculated fatigue lives for the fatigue GP-based model for the first training dataset and the following covariance functions: (**a**) Matern 3/2, (**b**) Matern 5/2, (**c**) RQ/SE.

It is concluded that the fatigue GP-based model trained on the dataset without the non-proportional loading paths cannot predict the material behavior of CuZn37 brass. The non-proportional loading induced an additional fatigue phenomenon in the CuZn37 brass compared to the phenomenon induced under proportional loading.

**Figure 7.** Comparison of error indexes obtained by application of different kernels to the fatigue GP-based model trained on the first training dataset.

The results obtained for the fatigue GP-based model trained on the second dataset, including one non-proportional loading path NPR, are displayed in Figure 8. By adding the results for one non-proportional loading path with 10 specimens to the training data, the prediction performance of the GP model was improved. An improvement is observed for all the applied covariance functions. The error indexes for the test dataset and analyzed kernels are practically equal (Figure 9) with *MPE* = [−0.2, −0.1]% and *T*(0.95) = [2.4, 2.5].

The training dataset with 30 specimens successfully trained the GP model for fatigue life prediction of CuZn37 brass subjected to proportional and non-proportional loading paths with test data of 50 specimens. One standard non-proportional loading path with a 90◦ phase shift provided sufficient information on the fatigue behavior of the GP model. Consequently, the fatigue life was successfully predicted for the four tested nonproportional loading paths (NPR1, NPR2, NPR3, NPR4) of different and complex shapes (Figure 2) and for one proportional loading path (Prop). Furthermore, the common complex measures of the non-proportional loading effect on fatigue life, such as non-proportional factors [79–82], integral approaches [66,72,82] and enclosing surface methods [83,84], are not needed.

**Figure 8.** *Cont*.

**Figure 9.** Comparison of error indexes obtained by the application of different kernels to the fatigue GP-based model trained using the second dataset.

### **5. Summary and Conclusions**

The Gaussian process (GP) was applied to build the fatigue model for the life prediction of CuZn37 brass subjected to various multiaxial loading paths. The physics-based input quantities in the form of stress and strain components on the critical plane of maximum shear served as predictors in the GP-based model. These quantities are invariant under the rotation of the coordinate system, and thus, the trained fatigue GP-based model is consistent with the invariance principles. Five stationary covariance functions for the GP with length scales assigned to each predictor (anisotropic kernels) were implemented and analyzed. Two Matern class kernels (M3/2, M5/2), SE, and RQ were preliminarily accepted for final validation based on the relevance analysis of predictors. The predictive performance of the GP-based model was verified using two training datasets. The first training dataset included only uniaxial and pure torsion loading paths, and the second dataset included one non-proportional loading path. The detailed conclusions are as follows.

• The relevance analysis of the applied input quantities for the fatigue GP-based model revealed that the maximum shear strain and normal stress on the plane of maximum shear are the most decisive factors for the life prediction of CuZn37 brass.


The GP-based model can effectively substitute parametric fatigue life prediction models if physics-based predictors consistent with invariant principles are applied. Further validation for different types of materials under the mean stress effect and a wider fatigue life regime is needed.

**Author Contributions:** A.K.: conceptualization, methodology, software, data curation, investigation, writing—original draft, writing—review and editing, visualization. D.S.: investigation, resources, writing—review and editing. Ł.P.: investigation, resources, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science Centre, Poland (Decision No. 2021/41/ B/ST8/00257).

**Data Availability Statement:** The data that support the findings of this study are openly available in the Mendeley data repository at http://doi.org/10.17632/7fbkf6y7gv.1 (accessed on 25 October 2021).

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