**1. Introduction**

In the production and operation of machines, especially the means of transport, a constant effort to reduce costs and energy consumption is employed. This goal is achieved, inter alia, by weight reduction, which can be achieved by lowering safety factors, topology optimization, or by replacing traditional materials with alternative, less dense ones. Moreover, the design process departed from the infinite-life design strategy, which requires the stress to be lower than the fatigue limit. Increasingly, the machines are designed according to the safe-life design strategy for the durability estimated by the designer. When the repair is very expensive (e.g., jet engines), a damage-tolerant design strategy is applied. Here, the service of a machine with diagnosed damage was accepted [1]. Therefore, the design of machines and structures in accordance with the savings trends and the recommendations of the latest strategies requires accurate fatigue life estimation methods accompanied by an uncertainty estimation of predictions.

The drive to reduce costs also applies to the design process. For example, the use of computational models should not require expensive experimental studies to obtain additional and more reliable material data. The high expenditure especially regards the design against material fatigue because the fatigue damage process is difficult to model. First, it consists of many fatigue crack development periods, such as crack nucleation, small crack growth, and macroscopic crack growth, and these periods differ in physical damage mechanisms also influenced by size and notch effects [2]. Each period is affected by different driving forces and requires different approaches, such as the crack nucleation and crack growth approaches [3]. Second, the fatigue behavior of machines and structures generally exhibits spectacular stochastic behavior [4–6]. This results from the variability of fatigue loadings, parts geometry, material properties, and microstructures. Such a significant level of difficulty in fatigue assessment has led to the development of new fatigue models, despite the numerous models that have been developed in the past few decades [7–10]. These models simplified the complexity of material fatigue by using empirical or semiempirical approaches to relate the primary loading quantities and material properties with the fatigue life. The empirical models adjust the parameters of the regression equation to the experimental data. Semi-empirical models combine fundamental physical principles and an empirical approach. Despite the use of some physical foundations to formulate models, the selected stress/strain quantities and material parameters are related by arbitrarily formulated parametric functions to obtain the best fit for the results of the experiment.

It is concluded that in the design process of engineering structures, the selection of the fatigue model is the primary problem. The choice should theoretically depend on the material properties (e.g., brittle or ductile), the type of loading (uniaxial or multiaxial, deterministic or random, proportional or non-proportional, stress ratio), and the range of deformation (elastic or elastic-plastic). In practice, owing to the multiplicity and complexity of models, the choice is most often limited by the level of the designer's knowledge. Conversely, owing to the cost and time-consuming nature of fatigue tests, the choice also results from the limited availability of data on the material and loading.

Based on the briefly presented problems of fatigue life prediction and existing solutions, a machine learning (ML) approach is proposed as a substitute for the semi-empirical fatigue models. The main advantage of this model is that it does not require the selection of a parametric (predefined) form of the fatigue model. The fatigue ML-based model should self-accommodate existing data and correctly reflect the fatigue behavior for testing data. Among several ML approaches, the neural network (NN) is widely used [11–13] in the fatigue field. According to Chen and Liu [12], NN modeling requires a large number of data, and there is no standard procedure for selecting an optimal NN architecture. One of the alternative ML approaches is the Gaussian process (GP) for regression with unique features that favors its application in fatigue life prediction of materials.


Existing studies on the use of GP to predict the durability of engineering machines and structures can be classified into approaches based on pure correlation of measured signals with progressive degradation of mechanical systems and hybrid approaches, including physics-based quantities.

In the first approach, the tested signals are not related to any model of the damage mechanism. Therefore, the results of this approach have mainly practical and predictive but not explanatory purposes. For example, Huchet et al. [21] used environmental parameters, such as wind speed and direction, in the fatigue assessment of wind turbine structures. In the research of Aye et al. [22] and Hong et al. [23], the prediction of the remaining life of bearings was directly based on the acquired vibration signals. Mohanty et al. [24,25] applied online signals from piezoelectric sensors attached to selected aircraft components in the GP for fatigue crack growth prediction. Hirvoas et al. [26] applied the GP to reduce uncertainties in a wind turbine numerical model considering different input data as support and blade structural properties.

Hybrid data-driven models involve input quantities related to the recognition of the failure mechanisms of a system. Sło ´nski [27] used the GP to identify concrete properties, among others, using minimal and maximal uniaxial stress values. Hu et al. [28] applied GP to estimate the uncertainty of fatigue crack growth in turbine discs in a study on fatigue crack growth evaluation. Ling and Mahadevan [29] proposed replacing computationally expensive finite element analysis for fatigue crack growth with the GP model. Farid [30] predicted fatigue failure under stochastic loading using a stress signal at the critical section of the mechanical component. The GP model was combined with an artificial neural network to enhance the predictive performance and provide uncertainty quantification.

In the briefly reviewed papers, the GP models were adjusted and trained for a given mechanical system. Consequently, the trained models cannot be applied to other mechanical systems, and thus, they can be mostly classified as health monitoring systems [31].

Karolczuk and Sło ´nski [32] proposed a novel approach in which the GP operates as a multiaxial fatigue life prediction model. The normal and shear stress amplitudes on the critical plane were selected as the physics-based input quantities for the GP model with the application of the squared exponential covariance function. The model was successfully verified on S355N steel and 2124 T851 aluminum alloy under the cyclic proportional combination of bending and torsion loadings at the high cyclic fatigue regime (stress-based condition). The proposed novel approach for fatigue life prediction requires further research and validation, especially under multiaxial non-proportional loading and stress–strain conditions.

This research aims to validate the GP applied to build a multiaxial fatigue model for the life prediction of CuZn37 brass under proportional and non-proportional loadings and stress–strain conditions.

The scope of this research involved an overview of selected classic (parametric) fatigue life prediction models (Section 2), basic concepts (Appendix A), and covariance functions (Appendix B) of the GP. Next, the fatigue life was estimated with GP using different covariance functions. As a result of the calculations, the physical quantities of key importance for the fatigue process were selected and compared with the quantities indicated in the analysis of classic models. The GP results were validated by comparing the results with the results obtained using the parametric models proposed by Fatemi-Socie, Brown-Miller, Glinka et al., and Yu et al. The calculations were conducted for eight loading cases, that is, axial, torsional, combined proportional axial-torsion loading, and non-proportional loadings. The level of non-proportionality was different in the case of the applied loadings owing to the phase shift (90◦ out-of-phase) and various frequencies of the applied strains (four different asynchronous loadings).

### **2. Brief Review of Fatigue Life Prediction Models**

Commonly applied semi-empirical fatigue models consist of fatigue damage parameters, which are scalar functions of spatial stress/strain components and reference fatigue curves. The fatigue damage parameter is used to reduce the spatial stress or strain state to a scalar quantity of the dimensions of stress, strain, or energy. It is then compared with the reference regression curves to calculate the fatigue life [33,34]. Uniaxial regression curves show an explicit relationship between the applied stress, strain, or energy values and the number of cycles to failure. These curves are the result of fatigue tests performed on a limited set of specimens and require the adoption of statistical assumptions [17]. Uncertainties in material parameters, loading and geometry can also be included in life prediction by the application of probabilistic modeling with sampling techniques [35–37].

The fatigue models differ in physical quantities, which were adopted as decisive factors for the fatigue process—predictors for the GP. For example, stresses, strains, loading non-proportionality factors, or strain energy, including elastic and plastic parts. These quantities, which fluctuate with time, are mostly reduced before incorporating them into the damage model to statistical parameters such as amplitudes, mean values, or maxima.

In this section, several selected models representing different physical principles are briefly described. A broader overview of the multiaxial fatigue models can be found in the papers concerning the critical plane approach [7], energy approach [38], or nonproportionality of loading [39].

### *2.1. Empirical Models*

An example of a purely empirical model can be the "ellipse quadrant" proposed by Gough and Pollard for ductile materials [40,41]. Because this is an equation for a particular loading case, that is, torsion and bending, a model is a function of the applied shear *τa* and normal *σa* stress amplitudes, as follows:

$$
\left(\frac{\tau\_a}{t\_{-1}}\right)^2 + \left(\frac{\sigma\_a}{b\_{-1}}\right)^2 = 1,\tag{1}
$$

where *t*−<sup>1</sup> and *b*−<sup>1</sup> are the fatigue limits for fully reversed torsion and bending, respectively. For brittle materials, Gough proposed the "ellipse arc" equation as follows:

$$
\left(\frac{\pi\_a}{t\_{-1}}\right)^2 + \left(\frac{b\_{-1}}{t\_{-1}} - 1\right)\left(\frac{\sigma\_a}{b\_{-1}}\right)^2 - \left(2 - \frac{b\_{-1}}{t\_{-1}}\right)\frac{\sigma\_a}{b\_{-1}} = 1. \tag{2}
$$

The above equations define the fatigue limit state of the applied stress amplitudes. It can be developed for a state at an arbitrary number of cycles to failure, as proposed in [42]. Empirical models have limited application (only combined torsion and bending loading) because they are not consistent with invariant principles.

### *2.2. Stress Invariants Models*

The two examples presented here are attempts to adopt the Huber–Mises yield criterion for fatigue by considering the hydrostatic stress. Both are a linear combination of the amplitude of the second deviator invariant *J*2, and the mean or maximum value of the hydrostatic stress *σH*, as in the case of the Sines [43] or Crossland [44] model, respectively:

$$
\sqrt{f\_{2,a}} + k\_{\text{s}} \cdot \sigma\_{H,\text{m}} = f\left(\mathcal{N}\_f\right), \\
\sqrt{f\_{2,a}} + k\_{\text{c}} \cdot \sigma\_{H,\text{max}} = f\left(\mathcal{N}\_f\right), \tag{3}
$$

where *ks*, *kc* are material parameters, *f Nf* is the reference fatigue curve.

The fatigue models based on stress invariants are criticized mainly because of problems in their implementation of random loading and non-proportional loading [7,45]. To overcome these problems, special procedures must be implemented.

### *2.3. Critical Plane Models*

Critical plane models are based on the observation of fatigue crack formation. Based on the observation, the fatigue cracks in metallic materials nucleate and develop in certain preferred planes within the material [46,47]. Thus, the critical plane approach assumes that the stress/strain components on a specific plane are primary for fatigue crack initiation and failure.

The stress-based Findley criterion is one of the first critical plane multiaxial fatigue models [48]. This is a linear combination of the maximum normal stress *<sup>σ</sup>n*,*max* and shear stress amplitude *<sup>τ</sup>ns*,*<sup>a</sup>* on the plane for which the equation reaches its maximum, as follows:

$$\max(\pi\_{ns,a} + k \cdot \sigma\_{n,\max}) = f\left(N\_f\right),\tag{4}$$

where *k* is a material parameter.

Brown and Miller [49] proposed a strain-based fatigue damage parameter composed of shear *γns*,*<sup>a</sup>* and normal *<sup>ε</sup>n*,*<sup>a</sup>* strain amplitudes on the plane experiencing the maximum shear strain amplitude. Kandil et al. [50] proposed a fatigue life prediction model based on this concept of the fatigue damage parameter, in the form

$$
\gamma\_{ns,a} + k\_{BM} \cdot \varepsilon\_{n,a} = f\left(N\_f\right),
\tag{5}
$$

where *kBM* is a material parameter.

Fatemi and Socie [51] proposed to relate the shear strain amplitude *γns*,*<sup>a</sup>* and maximum normal stress *<sup>σ</sup>n*,*max* normalized by the yield strength *<sup>σ</sup>yield*, in the following form:

$$\gamma\_{ns,a}\left(1+k\cdot\frac{\sigma\_{n,max}}{\sigma\_{yield}}\right)=f\left(N\_f\right).\tag{6}$$

Carpinteri et al. [52] proposed a nonlinear function of amplitude and mean value of normal stress and shear stress amplitude on the critical plane related to the average principal stress directions in the following form:

$$\sqrt{(\sigma\_{\rm n,a} + a\_{\rm CS}\sigma\_{\rm n,m})^2 + b\_{\rm CS}\tau\_{\rm ns,a}^2} = f\left(N\_f\right),\tag{7}$$

where *aCS*, *bCS* are material parameters.

Papuga–R ˚užiˇcka [53] also proposed a nonlinear function of amplitude and mean value of normal stress and shear stress amplitude on the critical plane of its maximum, as follows

$$\max\_{n} \left\{ \sqrt{a\_{PR} \tau\_{ns,a}^2 + b\_{PR} (\sigma\_{n,a} + c\_{PR} \sigma\_{n,m})} \right\} = f\left(N\_f\right),\tag{8}$$

where *aPR*, *bPR*, *cPR* are material parameters. The proposed formula was developed [45,54] to take into account the mean shear stress effect.

Glinka, Shen, and Plumtree [55] proposed the strain energy parameter accounting for both strains and stresses in the plane of maximum shear strain, as follows:

$$
\gamma\_{\rm ns,a}\sigma\_{\rm ns,a} + \varepsilon\_{\rm n,a}\sigma\_{\rm n,a} = f\left(N\_f\right). \tag{9}
$$

This model was modified by Pan-Chun-Chen [56] by introducing a weighting factor *kG* to normal components, as follows:

$$
\gamma\_{ns,a}\tau\_{ns,a} + k\_G\varepsilon\_{n,a}\sigma\_{n,a} = f\left(N\_f\right). \tag{10}
$$

Ince and Glinka [57] introduced a generalized strain energy fatigue damage parameter as a function of the elastic and plastic strain energy density contributed by the normal and shear stresses and strains on the critical plane of its maximum in the following form:

$$\max\_{\mathbf{n}} \left( \pi\_{\text{Its},\mathbf{a}} \gamma\_{\text{ns},\mathbf{a}}^{\varepsilon} + \pi\_{\text{Its},\mathbf{a}} \gamma\_{\text{ns},\mathbf{a}}^{\mathcal{P}} + \sigma\_{\text{n,max}} \varepsilon\_{\text{n},\mathbf{a}}^{\varepsilon} + \sigma\_{\text{n,max}} \varepsilon\_{\text{n},\mathbf{a}}^{\mathcal{P}} \right) = f\left( \mathcal{N}\_{f} \right),\tag{11}$$

where *n* is a unit vector that determines the orientation of the plane.

Yu et al. [58] modified the Ince-Glinka parameter and proposed the following model:

$$\frac{\gamma\_{\text{ns},a}\tau\_{\text{ns},a}}{\tau\_f'} + \frac{2\varepsilon\_{\text{n},a}\sigma\_{\text{n},\text{max}}}{\sigma\_{\text{yield}} + \sigma\_f'} = f\left(\mathbf{N}\_f\right),\tag{12}$$

computed on the plane of maximum shear strain, where *τf* and *σf* are material parameters deduced from the uniaxial reference curves.

### *2.4. Summary and Model Selection*

It is assumed that the GP-based fatigue model (Appendix A) can substitute any parametric functions proposed for fatigue life prediction. The substitution should be effective because the GP can map any fatigue behavior deduced from training data, and thus, the selection problem of adequate parametric function, for example, Equations (1)–(12), are omitted. The effectiveness of the GP-based model is the highest if the input data vector for the GP includes fatigue damage-related quantities. Based on the effort of many researchers and their attempt to fit the experimental data to parametric models, the primary quantities (predictors for the GP) for fatigue damage of metallic materials can be selected from the

shortlist presented in Table 1. The five predictors were selected from the critical plane models as being consistent with the observed fatigue damage mechanism of metallic materials. The role of these quantities and their interaction are discussed in [59,60]. The quantities from stress-invariant models could also be considered; however, owing to the ambiguous definition of amplitudes for stress invariants [61], they were neglected in the present research. Stress-based quantities are commonly applied at high cyclic fatigue (HCF), whereas strain-based quantities are common at low cyclic fatigue (LCF). Under non-proportional loading, the principal stresses rotate, which activates the higher number of slip systems and induces the complex interaction between dislocation movements. However, applying both quantities, i.e., strain and stress, the above mechanisms are reflected in predictors and the machine learning model is able to recognize this pattern to more effectively predict the fatigue life. The application of stress-and strain-based quantities reflects any possible additional material hardening effects occurring under non-proportional loading for some metallic materials. In the present research, we consider the plane or maximum shear strain amplitude appropriate for a wide class of metallic materials with dominant shear or mixed shear/tensile damage mechanisms.


**Table 1.** Primary predictors for fatigue life of metallic materials.

The fatigue life prediction performance of the proposed GP-based model was compared with the performance of four parametric fatigue models, including the proposal of Brown–Miller (5), Fatemi–Socie (6), Glinka et al. (8), and Yu et al. (10). The selected models required a more detailed description of the implemented reference fatigue curves *f Nf* and weighting material parameters. The basic regression curve implemented in defining the reference curve for each model is based on the Manson–Coffin curve [62] under fully reversed cyclic torsion loading, as follows:

$$
\gamma\_f \left( N\_f \right) = \frac{\pi\_f'}{G} \left( 2N\_f \right)^{h\_0} + \gamma\_f' \left( 2N\_f \right)^{c\_0}, \tag{13}
$$

where *τ f* is the shear fatigue strength coefficient, *G* is the shear modulus, *b*0 is the shear fatigue strength exponent, *γ f* is the shear fatigue ductility coefficient, and *c*0 is the shear fatigue ductility exponent. The Brown–Miller, Fatemi–Socie, and Glinka et al. models include the weighting material parameters that can have a life-dependent form [33,63] to be fully consistent with two uniaxial strain-life curves obtained under tension-compression and torsion loadings. The torsion strain-life curve is given by Equation (13), and the tension-compression strain-life curve is

$$
\varepsilon\_f \left( N\_f \right) = \frac{\sigma\_f^{\prime}}{E} \left( 2N\_f \right)^b + \varepsilon\_f^{\prime} \left( 2N\_f \right)^c \tag{14}
$$

where *σf* is the axial fatigue strength coefficient, *E* is Young's modulus, *b* is the axial fatigue strength exponent, *εf* is the axial fatigue ductility coefficient, and *c* is the axial fatigue ductility exponent. The first component of Equations (13) and (14) are the elastic strain and plastic strain, respectively. Decomposition into elastic and plastic parts is consistent with the physical mechanism, and it is necessary to derive life-dependent material weighting factors appropriately. This decomposition leads to the following notations:

$$\varepsilon\_f^\varepsilon \left(\mathbf{N}\_f\right) = \frac{\sigma\_f^\prime}{\mathcal{E}} \left(2\mathbf{N}\_f\right)^b,\\ \varepsilon\_f^p \left(\mathbf{N}\_f\right) = \varepsilon\_f^\prime \left(2\mathbf{N}\_f\right)^c,\\ \sigma\_f \left(\mathbf{N}\_f\right) = \sigma\_f^\prime \left(2\mathbf{N}\_f\right)^b,\\ \tau\_f \left(\mathbf{N}\_f\right) = \tau\_f^\prime \left(2\mathbf{N}\_f\right)^{b\_0},\tag{15}$$

where the upper indexes *p* and *e* indicate the plastic and elastic strain parts, respectively. Implementing these notions, the detailed formulas for the parametric life prediction models were obtained [63], as follows:

•Brown–Miller model:

$$f\begin{pmatrix} N\_f \\ \end{pmatrix} = \gamma\_f \begin{pmatrix} N\_f \\ \end{pmatrix}\_{\text{s.t.}} \tag{16a}$$

$$k\_{BM}\left(N\_f\right) = \frac{2\left[\gamma\_f\left(N\_f\right) - (1+\nu^c)\varepsilon\_f^c\left(N\_f\right) - (1+\nu^p)\varepsilon\_f^p\left(N\_f\right)\right]}{(1-\nu^c)\varepsilon\_f^c\left(N\_f\right) + (1-\nu^p)\varepsilon\_f^p\left(N\_f\right)},\tag{16b}$$

where *ν<sup>e</sup>* and *νp* are the elastic and plastic Poisson ratios, respectively.

•Fatemi–Socie model:

$$f\left(N\_f\right) = \gamma\_f\left(N\_f\right) \tag{17a}$$

$$k\_{FS}\left(N\_f\right) = \left[\frac{\gamma\_f\left(N\_f\right)}{(1+\nu^e)\varepsilon\_f^e\left(N\_f\right) + (1+\nu^p)\varepsilon\_f^p\left(N\_f\right)} - 1\right] \frac{2\sigma\_{yield}}{\sigma\_f\left(N\_f\right)}\tag{17b}$$

• Glinka et al. model:

$$f\left(N\_f\right) = \gamma\_f\left(N\_f\right)\tau\_f\left(N\_f\right) \tag{18a}$$

$$k\_{\mathbb{G}}\left(N\_{f}\right) = \frac{4\gamma\_{f}\frac{\pi\_{f}\left(N\_{f}\right)}{\sigma\_{f}\left(N\_{f}\right)} - 2\left((1+\nu^{\varepsilon})\varepsilon\_{f}^{\varepsilon}\left(N\_{f}\right) + (1+\nu^{p})\varepsilon\_{f}^{p}\left(N\_{f}\right)\right)}{(1-\nu^{\varepsilon})\varepsilon\_{f}^{\varepsilon}\left(N\_{f}\right) + (1-\nu^{p})\varepsilon\_{f}^{p}\left(N\_{f}\right)}\tag{18b}$$

• Yu et al. model:

$$f\left(N\_f\right) = \gamma\_f' \left(2N\_f\right)^{c\_0} + \frac{\mathbf{r}\_f'}{G} \left(2N\_f\right)^{2h\_0}.\tag{19}$$

Substituting Equations (13)–(19) in (6) and (8)–(10) and implementing an iterative gradient-based procedure results in the calculation of the fatigue life *Ncal* for each parametric model.
