An Approach to Predicting Fatigue Crack Growth Under Mixed-Mode Loading Based on Improved Gaussian Process

This paper proposes an approach to predicting fatigue crack growth under mixed-mode loading based on improved Gaussian process. In terms of analyzing the theoretical background for fatigue crack growth, a corresponding finite element model is built to generate sufficient simulation data, which is utilized to obtain the key parameters (e.g., stress intensity factor) for the fatigue crack growth process. And then, a Gaussian process model is achieved to meet the condition that the stress intensity factor is a nonlinear continuous change in crack growth, especially for mixed-mode loading. Following, an idea of local sample densification method is implemented to improve Gaussian sample generation process according to the simplified model of the crack growth path. Based on the above investigation, a fatigue crack growth prediction model using the improved Gaussian process is finished, which is subsequently verified through the test data of lower bainite steel (SCM435) material. The results show that the proposed approach has better computational accuracy and efficiency than the traditional finite element method in predicting fatigue crack growth under mixed-mode loading.


I. INTRODUCTION
It is well known that engineering structures are prone to fatigue crack under cyclic stresses during working pro-cesses. Once fatigue crack accumulates to a certain level, it may cause safety accidents of structures. Taking the collapse of the Kielland rig at the Ekofisk field in the UK North Sea in 1980 as an example, its fatigue cracks produced in individual support legs continued to grow under loads in the marine environment, which caused the entire platform overturned eventually [1]. This accident led to significant economic losses and casualties. A reasonable approach to avoiding this kind of accident is to predict fatigue crack growth for determining some precautions. Therefore, fatigue crack growth prediction has attracted many scholars' attention because The associate editor coordinating the review of this manuscript and approving it for publication was Agustin Leobardo Herrera-May . of the meaningful engineering application and scientific significance.
The general method of analyzing and predicting fatigue crack growth for engineering structures is using the theory of fracture mechanics, which indicates crack growth through clarifying the corresponding mechanism. To investigate the inherent pattern of crack growth and obtain the approach to calculating crack growth rate, many researchers have performed relevant theoretical and experimental studies. For example, Irwin [2] considered the singularity of crack tip stress and proposed concept of stress intensity factor (SIF), which was a parameter related to the load size, load angle, crack geometry and crack length. Paris [3] gave an exponential relationship expression based on SIF to calculate crack growth rate. This is the so-called Paris equation, which is well-known and widely accepted by researchers. But this equation only explains rules of crack growth rate in stable region. To overcome the drawbacks, Forman et al. [4] VOLUME 9,2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ added the fracture toughness to Paris equation, thus this equation can be applied to the region with high-speed growth. Donahue et al. [5] proposed a threshold value of SIF range to deal with the nonlinearity profile for the crack growth rate curve. To further consider the crack closure phenomenon, Elber [6] improved this equation to analyze direction and rate of crack growth. But this equation still should subject to constant amplitude loading (CAL) condition. When fatigue crack is under the condition of variable amplitude load (VAL), its growth rate shows a great change. This phenomenon is called overload retardation phenomenon. Lu et al. [7] conducted a series of experiments on automobile steel under different load ratios and stress ratios, and then proposed a new SIF calculation method. Ding et al. [8] carried out single overload and underload tests on the central crack specimen of Q345R steel. The results showed that crack growth rate would be reduced or increased in a short time with overloaded or underloaded conditions. To describe this phenomenon, many researchers have proposed different theoretical explanations and improved formulas. For example, Wheeler [9] proposed that the overload retardation occurred because overload created a large plastic region at the crack tip. Currently, a crack growth rate model proposed by McEvily et al. [10] is widely accepted. This model can not only explain the phenomenon of overload retardation but also stimulate the effect of initial defects. But this modified model is only applied in linear and near threshold regions. Cui et al. [11], [12] improved McEvily's model to make it suitable for the condition of instability fracture. Although there have been many achievements that can be utilized to analyze and predict fatigue crack growth, the proposed formulas in the above literature mainly aim at crack growth rules under mode-I loading. Because of the structure in general subjects to mixed-mode loading in practical operation, it is insufficient to only investigate fatigue crack growth with mode-I. Under the influence of this fact, the crack growth path will produce deflection, which greatly increases difficulty of calculating SIF at crack tip. Focusing on this situation, Richard et al. [13] designed a composite load loading device to perform fatigue tests on compact tension shear (CTS) specimen, and then gave an analytical solution of SIF of the I-II composite crack tip at different loading angles. But this formula requires the crack length to be within the specified range. Thus, it does not obtain a SIF for deflected crack tips. Therefore, the structure cannot be predicted by mathematical modeling methods for fatigue crack growth under mixed-mode loading. Besides, some researchers calculated the SIF at the tip of a deflected crack by direct observation approaches. Silva et al. [14] used digital image correlation (DIC) to directly extract the displacement field of the crack geometry and calculate the SIF through image processing. Although the results showed that this method gave more accurate results, it was too complicated and expensive to be used in practical applications.
Nowadays, due to the complexity of crack growth under mixed-mode loading, there is not a unified analytical prediction model of fatigue crack growth. With the rapid development of computer technology, the finite element method (FEM) has shown great applicability and flexibility. Therefore, the FEM-based fatigue crack growth prediction method has been gradually proved to an effective numerical solution. Maligno et al. [15] developed a 3D numerical analysis method for shaft structures based on the FEM. This method can solve the crack tip SIF calculation problem by continuously dividing meshes according to the set of crack growth increments a. Fixed-SIF and Linear-SIF are two specific application methods for calculating the crack growth rate. The former considers the SIF calculated at each growth step as a fixed value, while the latter brings the SIF calculated from the current growth step and the next growth step to the linear fitting formula into the calculation process. But Fixed-SIF has the problem of discontinuity of crack growth in the calculation process, and it can only bring the finite element results of the SIF in the current crack state into the mathematical model for calculation. Thus, this method has the disadvantage of low calculation accuracy. Regarding Linear-SIF, it assumes a linear relationship between the SIF of different crack lengths, which is inconsistent with the actual situation. The traditional FEM requires constant re-meshing, and this increases workload and reduces computational efficiency. The extended finite element method (XFEM) has the feature of avoiding re-meshing for fatigue crack growth simulating processes. Thus, it solved the calculation discontinuity problem of FEM. Pandey et al. [16] simulated the law of high-cycle fatigue crack growth by combining the continuous damage mechanics method with XFEM. Xi et al. [17] used XFEM to simulate the crack growth of rock and discussed the influence of crack surface friction on crack pattern and growth. But the current XFEM is based on the principle of linear elastic fracture mechanics (LEFM), the influence of the delay effect on the crack growth cannot be considered. Dirik et al. [18] used XFEM to consider the problem of crack growth under variable amplitude loads. The process is mainly to calculate crack growth rate by using the results of the finite element calculation of the SIF into the correction formula. Due to the superiority of XFEM in terms of mesh reconstruction, this method improves the computational efficiency. However, the prediction accuracy is still limited to a set amount of crack growth.
In recent years, with the development of machine learning technology, many researchers have proposed to use a data-driven approach to solve the crack growth prediction problem. Since the Gaussian Process (GP) regression is an effective tool that only uses simple linear algebra to deal with nonlinear models, it is often used in such studies. Neerukatti et al. [19] combined a data-driven model with a physical model to achieve an accurate prediction of mode-I crack growth results. And in their subsequent study, GP was used to establish the mapping relationship between the crack tip energy release rate G and the crack extension rate to improve the prediction accuracy under biaxial loading [20]. Yuan et al. [21], [22] used GP to establish the mapping relationship between crack length and health monitoring signal and corrected the errors arising from the prediction process. However, there is a lack of research on data-driven fatigue crack growth prediction methods for general structures under mixed-mode loading. Because SIF is related to crack growth paths and the structure may have numerous crack growth paths under unknown boundary conditions. Thus, to obtain accurate fatigue crack growth prediction results, the SIF-GP model will have a large amount of sample data, which also leads to long model calculation time and large fitting errors. There are two ways to solve this problem: 1) establish sparse sample space, as Cai et al. [23] used a sparse Gaussian model for the computer vision recognition of human activity and Kuzin et al. [24] used this method for online Bayesian inference; 2) establish local GP model, as Gramacy et al. [25] proposed method of divided the data to establish multiple local GP models in place of the global GP model. The former has the problem of prediction accuracy, while the latter needs to study how to divide the data and allocate it to different models. Thus, how to improve the accuracy of the GP model while maintaining a reasonable sample size becomes a new problem.
Therefore, to solve the problem of fatigue crack growth prediction described above, this paper proposes an approach for predicting fatigue crack growth based on the GP model with local sample densification. Using FEM to get the sample space of SIF, and then a proxy model of the SIF is established based on the GP. To control the number of samples to improve model computational efficiency, this paper builds an improved generation approach of sample space. First, a sparse sample space is generated in the possible crack growth area in the structure. Subsequently, the sample is densified according to the real-time prediction result or the actual crack growth situation, which not only ensures the reasonableness of the sample number but also improves the accuracy of the prediction. The SIF-GP model and the fatigue crack growth model are combined to calculate the crack growth rate, and the deflection angle of the crack growth is calculated based on the maximum circumferential stress criterion. Finally, this paper presents the prediction of experimental data of lower bainite steel (SCM435) in the literature and compares the accuracy of the prediction results with the traditional FEM. It is proved that the approach proposed in this paper can greatly improve the efficiency and accuracy of fatigue crack growth under mixed-mode loading.

II. THEORETICAL BACKGROUND
To establish a SIF-GP model, the authors need an approach that can accurately and quickly produce SIF samples. With the development of computer technology, numerical analysis methods based on FEM have shown great advantages. The characteristic of the method is that the calculation process is not restricted by geometric conditions. Therefore, based on the commercial software ANSYS, this section builds the finite element modeling approach of the CTS specimen and the calculation method of the SIF at the crack tip. Meanwhile, to predict the crack growth path and growth rate, the prediction criterion of crack growth deflection angle and the related theory of growth rate model are introduced subsequently.

A. CTS SPECIMEN MODELING METHOD
To study the fatigue crack growth prediction approach under mixed-mode loading, the CTS specimen and fixture are used as the research objects. Its geometric structure is shown in Fig. 1, and the specimen can be loaded from mode-I to mode-II. For example, the loading is pure mode-I when the angle is 0 • , the loading is pure mode-II when the angle is 90 • , and the loading is mixed-mode when the angle is between 0 • and 90 • . To simplify the calculation process for improving the calculation efficiency, this paper assumes that the research object is penetrating crack. Thus, the PLANE183 element was used for meshing and the 3D crack shape characteristics are not considered. Meanwhile, according to the influence of each part of the model on the SIF results, the mesh sizes are divided into three types, as shown in Fig. 2. Since the fixture is not the key part of the calculation, it is coarsely meshed for improving the efficiency of the calculation. To keep the accuracy of the calculation results, the specimen is divided into fine mesh, and singular element is set at the crack tip.
As shown in Fig. 3, the connection between the CTS specimen and the fixture is achieved by coupling the nodes where the bolt holes are located. Meanwhile, by establishing coupling constraints between the center of the loading hole and the circumferential surface, the load can be equivalently applied at the center of the circle.

B. INTERACTION INTEGRAL METHOD
After the finite element model is established, a method needs to be selected to obtain the SIF. The interaction integral method is widely used because it is independent of the path and can bypass the singular area of the crack tip. It is based on the J integral to separate and obtain SIF in the real field by establishing an auxiliary field at the crack tip.
Therefore, J integral defines the lower loop: where u i is the component of the displacement vector, W is the strain energy density, W = ε ij 0 σ ij dε ij , T i is the tension vector and T i = σ ij n i , is the arbitrary smooth closed loop, n i is the unit normal vector of the integral loop.
The relationship between J integral and SIF is as follows.
It is assumed that the component receives the real load and auxiliary load. Bring the superposition of the real field (σ a ij , ε a ij , u a ij ) and auxiliary field (σ a ij , ε a ij , u a ij ) under the load into (1), and the J integral obtained by the two fields alone are removed. The interaction integral between the two states can be obtained by the following expression: According to (2), the result can be equivalent to Bring K a I = 1, K a II = 0 and K a I = 0, K a II = 1 into the formula.
where I 1 and I 2 are the interaction integral value in two cases, respectively.

C. METHOD FOR DETERMINING CRACK GROWTH DEFLECTION ANGLE
The crack growth path has a direct effect on the SIF. When the structure is subjected to mixed-mode loading, the path will be changed. Currently, there are three main criteria for determining crack deflection angles: 1) Maximum Energy Release Rate (MERR) criterion proposed by Hussain and Pu [26]; 2) Maximum Tangential Stress (MTS) criterion proposed by Redogan and Sih [27]; 3) Maximum Strain Energy (SED) criterion proposed by Sih [28]. Based on Floros et al. [29] study of the predictive power of different crack growth direction criteria through finite element simulations, this paper selects MTS to simulate the crack deflection path, and the criterion assumes that the crack growth direction is determined by circumferential stress σ θ θ . It can be expressed as follows: As long as the circumferential stresses meet ∂σ θ θ ∂θ = 0 and ∂ 2 σ θθ ∂ 2 θ < 0, the first-order partial derivative solution can be obtained The crack deflection angle σ 0 can be written as: As can be seen from the above equation, K II determines the direction of deflection of the crack. The sign of the crack deflection angle θ 0 depends on the positive and negative signs of the K II , as Fig. 4. When K II is positive (shearing force is positive according to the elastic mechanics' sign, K II is positive.), θ 0 is negative (clockwise around the local coordinate system x .). Conversely, when K II is negative, θ 0 is positive (anticlockwise around the local coordinate system x .).

D. FATIGUE CRACK GROWTH RATE MODEL
A crack growth rate model that can describe the crack growth process is crucial to obtain accurate predictions of the crack growth life. In this study, a modified formula based on the Paris model proposed by Huang [30] will be used. This modified model can consider the influence of the stress ratio R on the fatigue crack growth rate, and the formula is expressed as follows: where C is a material and environmentally constant of dimensions, m is a crack growth law exponent, K th0 is a threshold SIF range at R = 0 (MPa √ m), M R is a correction factor for the loading ratio, M P is a correction factor for the loading sequence interaction.
In curved crack growth, the equivalent SIF K eq at the crack tip has to be used for accurate life estimation. K eq are suggested by Tanaka [31], and he used in the equation instead of K .
For the retardation phenomenon, the Wheeler model is usually used to describe its fatigue crack growth behavior. But it is found that the underload after the overload will weaken the retardation effect in practical applications. If the underload occurs first, there is no obvious effect on the crack growth rate. Therefore, the model introduced the concept of the underload plastic zone and revised the retardation coefficient in the Wheeler model, as shown in Fig. 5. The improved model can consider the influence of the load sequence, such as underload, overload-underload, and underload-overload. The correction factor for the loading sequence interaction M P is defined as follows: VOLUME 9, 2021 where n is an exponent used to accelerate or retard crack growth due to load interactions, α is a plastic zone correction factor, a OL is the length of the crack at overload, r OL is the size of the plastic zone of the crack tip at overload, r is the size of the plastic zone of the crack tip at underload, r y is the plastic zone radius for the current crack geometry, K maxOL is the maximum SIF at overload, K minCAL is the minimum SIF at constant load, K minUL is the minimum SIF at underload.

III. IMPROVED GAUSSIAN PROCESS MODEL
As described in Part I, in previous studies of crack growth under mixed loads, it is necessary to bring the crack information in the current state back to the finite element model to calculate the SIF. This method causes low computing efficiency and resource waste. GP regression is essentially a non-parametric kernel method based on a Bayesian framework that uses statistical principles to establish a mapping relationship between a dependent variable and an independent variable. Therefore, this part will first introduce the theory of GP regression. Subsequently, the improved SIF-GP model for mixed-load is then proposed based on a simplified model of the crack growth path.

A. GAUSSIAN PROCESS REGRESSION THEORY
Suppose that the training data set is Considering the existence of errors and noise in the observed data, the relationship between the observed output y and the true output f (x) are defined as: where ε is white Gaussian noise, and it is a normal distribution consistent with zero mean and unit variance. If f (x) constitutes a Gaussian random process, its mean function µ (x) and covariance function k x, x are defined as: For simplicity, only the zero mean µ (x) = 0 can be considered. If f (x 1 ), . . . , f (x n ) follows a Gaussian distribution, the joint distribution of noise observation vector y is also a zero-mean. It can be expressed as follows: where I n ∈ R n×n is the unit matrix, N µ, is the Gaussian distribution of the mean vector µ and the covariance matrix , n × n is the covariance matrix K(X, X), whose elements consist of k x i , x j , i, j = 1, 2, . . . , n.
The joint probability distribution formed by the target output vector y obtained from the subset of training data x and the corresponding output f * obtained from the test data x * can be expressed as: where k (X, x * ) = k T (x * , X) = [k (x 1 , x * ) , · · · , k (x n , x * )] T is the n × 1 dimensional covariance function vector between all training input data and the test input k (x * , x * ) is the selfcovariance function value of the test data. Therefore, the predicted mean and variance: GP model is determined by the prior covariance and the super-parameters. The prior covariance determines the properties of the GP model and plays a crucial role in the predictive effect of the model. To better capture the various characteristics of the data, this study uses an automatic correlation determination (ARD)-based covariance function to construct the GP model, as shown in Table 1. ARD can specify a separate length-scale l d for each input element, and l d reflects the relative importance of different input variables [33] (e.g., the variables with smaller length scales have a greater impact on the prediction). Thus, this feature makes the GP model greater predictive flexibility.

B. SIMPLIFIED MODEL OF CRACK GROWTH PATH UNDER MIXED LOAD
Generating suitable and accurate sample data is essential to build the SIF-GP model. When the structure is subjected to a mixed-mode loading, the crack growth path is deflected. However, as the crack continues to growth, the deflection angle of the path decreases gradually. Therefore, the crack growth path is not a straight line, which brings challenges to the definition of the input parameters of the SIF-GP model. As shown in Fig. 6, the loading direction is perpendicular to the x-axis. The initial crack l OA deflection angle is θ 2 under the load, and it grows a trace d a to point B. Meantime, when going through multiple growth steps, building a crack finite element model that follows the actual crack growth path exactly will greatly increase the difficulty and workload. Therefore, to predict the direction of the subsequent growth of the actual crack l OAB , this paper uses a simplified model of the crack growth path. This model equates the actual growth path l OAB to l OB . Thus, the equivalent crack length l OB and the crack deflection angle β can be respectively expressed as follows:

C. ESTABLISHING SAMPLE SPACE OF SIF
To obtain the SIF sample space, it is necessary to analyze the characteristics of the simplified model of crack growth paths under mixed-mode loadings developed above. Then, according to the path characteristics, the sample input format of the SIF-GP model is defined to complete the mapping relationship between the crack growth path and the SIF at the crack tip. In CTS specimens, when the load angle of the structure is changed, the crack growth angle will change significantly. Therefore, based on the simplified model of the crack growth path developed in the previous section, this paper assumes that the crack growth path at this loading angle is a straight line when the loading angle is not changed. In the case when the load angle of the structure changes multiple times during the entire loading cycle, the crack growth path also changes. The schematic representation of the crack growth path for this case is shown in Fig. 7.
Based on the number of load loading angle changes that occur during crack growth, this paper defines the SIF-GP sample format as follows: where A is the crack growth length increment under the current loading angle, n is the number of loading angle changes, i is the number of crack growth path information. In particular, when the structure has a mode-I crack growth at the beginning, the crack deflection angle α 0 is 0, and the sample space input is only the crack growth increment A 0 . Only if the model predicts or the actual crack detects a deflection in the crack path, then the sample space input changes to the form of (A n , α n ).

D. GP MODEL IMPROVEMENT APPROACH
To improve the prediction accuracy of the GP model, the previous approach increases the sample points to raise the sample density. With the continuous increase of sample points, this approach can consider all possible crack growth path conditions. However, the input of the SIF-GP model under mixed-mode crack is a high-dimensional vector, so the enumeration method is generally used to generate samples containing various variables to consider all possible situations.  Taking two angle changes in crack growth as an example, the sample point generation approach is shown in Fig. 8. First, according to the possible loading angle variation range and the preset maximum crack growth length, the variation range of each input parameter in the sample space is divided in equal proportions. The result is to output n 4 crack path information. Then, the information on the crack growth path represented by each sample is brought into the finite element software to simulate and calculate the corresponding SIF solution.
It is obvious that as the step size divided by each parameter becomes finer, the accuracy of the prediction results will increase. Meanwhile, when the load on the structure changes, this will add two variables to the sample space, which results in a huge increase in the number of samples. Although this approach can completely cover all possible crack growth path information, the number of samples obtained in the case of multiple angle changes is huge. The result not only leads to a huge cost burden due to a lot of early calculation sample data, but the excessive amount of data will make the unacceptable cubic complexity of the GP model. Therefore, the standard GP model needs to be improved to make it suitable for learning and prediction in the case of high-dimensional variables.
From the perspective of crack growth prediction, the time point of the change in the size and direction of the load on the structure is unknown, so there are many possibilities for the growth path. However, from the perspective of the entire crack growth cycle of a specific structure, due to the change in the load angle at a certain point during the growth period, the actual crack growth path is unique. When using the SIF-GP model to predict fatigue crack growth, only the sample points around the crack path obtained by prediction calculation or actual detection will contribute to improving the accuracy of SIF prediction. Most of the sample points at other locations have little impact on the prediction results, so it greatly wastes computing resources and increases computing time.
To solve this problem, this paper builds an improved GP model of SIF based on the local sample densification. The sample generation process of the improved model is shown in Fig. 9. The crack first grows when the loading angle is α 0 , such as C1. When the loading time is N1, the loading angle of the crack will change to α 1 , and the parameter A 0 in the sample space at this time is a fixed value. Subsequently, the deflection angle can be predicted by the maximum circumferential stress criterion. By adding a given range of deviation to this prediction, the possible path range under the current load angel is obtained, such as C2. The sparse sampling method is performed within this range, and the initial sample number is reduced by increasing the interval between sample points. According to the prediction or actual observation of the crack growth path, to improve the calculation accuracy of the SIF, random sampling and sample densification are performed around the path that may follow, such as C3. Similarly, when the loading time is N2, the loading angle changes again, and the parameters A 0 and (α 1 , A 1 ) in the sample space at this time also become certain values. Then repeat the above process to obtain the possible path range under the loading angle α 2 , and according to the predicted path conditions, the sample space is locally encrypted to improve the prediction accuracy of the stress intensity factor under the current crack path.

IV. FATIGUE CRACK GROWTH PREDICTION FRAMEWORK COMBINED WITH IMPROVED GP MODEL
To combine each of the above approaches, so this paper establishes the fatigue crack growth framework that can predict the crack growth path and growth life under the condition of mixed-mode loading. Firstly, the authors propose an improved GP-SIF model based on the local sample densification, which not only ensures the accuracy of the calculation results but also greatly reduces the calculation time of the sample space. Secondly, by combining the SIF-GP model with the mode-I crack growth rate model, the calculation efficiency of the crack growth results under mixed loading is improved. The overall program is shown in Fig. 10, and the steps of the program can be described in detail as follows: (1) A multi-scale FEM is developed for a specific structure, which contains different crack information, such as the initial length of the crack, the crack growth length, the crack growth angle, load size and direction. The geometric information and boundary conditions of the structure are regarded as the fixed quantity of the model. Use ANSYS parametric design language (APDL) [34] to write a program that can realize automatic modeling and analysis through modifying a few parameters. When the loading angle changes during crack growth, the program generates a sparse sample space based on the current crack growth, and then brings sample space into the finite element model to obtain the SIF. Combine the sample space with the kernel function to find the optimal super-parameters to build a SIF-GP model to predict the crack growth path range under the current load. To densify the samples to improve the prediction accuracy, based on the prediction results and actual crack observation results, randomly generate samples within the crack growth path area. By combining the SIF-GP model under the current loading angle with the crack growth rate model, the prediction of crack growth life and path under mixed-mode loading is completed. (2) According to the input load information, the program needs to analyze the variation of the load at the current moment. If there is overload or underload during the loading process, the retardation effect needs to be considered, and the d a/d N at the current time is corrected by using the correction factor for the loading sequence interaction M P . If the load keeps constant, the program will calculate the crack growth a i and the deflection angle θ i−1 . According to the crack path simplified approach, the initial crack length a i and the initial growth deflection angle θ i at the next moment are output to the program. The same procedure is followed in each increment, and when the pre-determined total length of crack growth is reached, the program will automatically output the coordinates of the crack growth path and predict crack growth life. It is noted that when the calculated K eq is greater than the fracture toughness of the material, the program will immediately stop the calculation and output the relevant result information. In the following section, we will use a specific case study to validate the effectiveness of the proposed approach in this paper.

A. DATA SOURCES AND LOADING CONDITIONS
Kim et al [35] experiments on fatigue crack growth under mixed-mode loading in lower bainite steel (SCM435) CTS specimens. With an excellent combination of strength and toughness, this steel is widely used in oil platforms, aerospace and other fields. The nominal composition of the steel is Fe-0.35C-0.21Si-0.78Mn-0.98Cr-0.18Mo in wt%. The thickness of the CTS specimen is 3.75 mm, the width is 45 mm, and the initial cut depth is 20 mm. The specific dimensions are shown in Fig. 11. The test frequency f = 10 HZ, and the load ratio R = 0.1.
Before the experiment, the CTS sample was preloaded with the crack in mode-I, so that the crack generated and grew forward 2 mm. The specimen is fixed to the lower fixture, while the upper fixture is connected to the fatigue testing machine. As shown in Fig. 1, by applying loads of  431.6N-4316.4N to the specimen in the directions of 0 • , 15 • , 30 • , 45 • and 60 • , data on the crack growth path and rate are obtained. This paper uses these data to compare and verify the fatigue crack growth prediction method based on improved GP model.

B. BUILDING IMPROVED GP MODEL
According to the SIF-GP model establishment approach, the first step is to obtain SIF samples at the crack tip under different crack paths. Since the load angle of each experiment does not change, the path information of the crack is only a group of crack length and crack deflection angle. Therefore, the input data format of the GP model is (β 1 , a 1 ). The construction procedure of the GP model is as follows: Step I: First, a sparse sample space is created based on the geometry of the structure. According to this specimen structure, the range of crack growth lengths is 22 mm-34 mm, and the range of crack deflection angle is 0 • -60 • . The sampling interval of crack length and deflection angle is 1 mm and 15 • , respectively. As a result, 60 sample points can be obtained to form GP sparse sample data, and these samples are uniformly distributed within the design range. Subsequently, according to the loading angle of each test, the crack growth deflection angle is predicted in terms of the MTS. Add an error value to the calculation result and set the value to 10 • to establish the possible area of the crack growth path. Finally, a total of 130 random sample densification sites are obtained in this region.
Step II: According to the finite element modeling approach described in the second section, the corresponding APDL program is achieved. By using the script language, each combination of crack information is brought into the finite element to calculate the corresponding K I and K II .
Step III: Sample results obtained by finite element calculation are divided into two parts. All sparse sample points and 100 densification samples are used to train the SIF-GP model, while the remaining 30 samples are used to determine the suitable covariance kernel function and verify the accuracy of the prediction results.
The five kernel functions of covariance listed in Table 1 are brought into the SIF-GP model for comparison. To quantitatively analyze the difference between the prediction results and find the covariance kernel function with the smallest error, this paper uses the root mean square error (RMSE) value as the evaluation index of prediction accuracy, and a smaller RMSE indicates better predictive power of the covariance kernel function [36].
where N is the number of samples, K i (x) and K i (x) are the i-th crack SIF finite element calculation value and GP prediction value, x is the type of SIF (mode-I or mode-II).  Table 2 shows the RMSE values of K I and K II when the loading angle is 30 • using the covariance kernel function based on ARD. The RMSE of K I and K II calculated using the ARD-Matérn (5) kernel function is the smallest, so ARD-Matérn (5) is selected as the covariance kernel function for calculating K I and K II . Besides, Fig. 12 shows a 3D graph of the fitted GP using the ARD-Matérn (5) kernel function and the testing samples computed from the FEM. It can be seen from the figure that the predicted values of GP are all close to the FEM values, and the maximum error when calculating K I does not exceed 1.5%, and the maximum error when K II does not exceed 1.2%. This result proves that the accuracy and reliability of the GP for predicting the SIF can well capture the nonlinear characteristics of the SIF under different crack information.

C. DISCUSSION ON PREDICTION RESULTS
The SIF calculation part in the fatigue crack rate model is replaced by the established SIF-GP model, which overcomes the problem that the original model cannot calculate the fatigue crack growth rate in the mixed-mode. Based on the experiment data, the crack growth rate da/dN = 8.1389 × 10 −12 K eq0 2.781 − 7.1 2.781 of the CTS specimen at a loading angle of 0 • was fitted. Fig. 13 shows that the result can fit the fatigue crack growth rate under various loading angles well. In Fig. 14, the curve of the experiment data of the crack growth at 0 • is compared with the crack growth prediction curve calculated by GP model prediction after fitting, where the N coordinate is the loading period and the a coordinate is the crack length. The maximum error is less than 0.2 mm, which shows that the crack growth rate model predictions have an excellent fit with the experiment results.
To verify the effectiveness of the MTS for the initial crack deflection angle, the experimentally measured initial crack   deflection angles under different load angles are compared with the deflection angles calculated by the MTS, as shown in Table 3. The error of the crack deflection angle is within the acceptable range for all the other loading angles except for 12.5% error in the crack deflection angle at 45 • . This result indicates that the MTS is applicable in predicting the mixed-mode loading case.  Fig. 15 is a graph showing the N-a curve calculated by different methods at each loading angle and the trend plot of the relative error between predicted and experiment values with increasing real experiment crack length. Fig. 16 shows the comparison of the calculation time of different methods. As can be seen in the figure, the accuracy of the crack growth data predicted by the Fixed-SIF with a = 0.6 mm is the lowest. As the number of loading increases, the error from the experimental value becomes larger and larger. At each loading angle, the final error will reach 4.0 mm-5.5 mm. The Fixed-SIF with a = 0.3 mm has a smaller crack growth increment, so its calculation accuracy is improved by about 20% compared the Fixed-SIF with a = 0.6 mm, but the overall calculation time is 100% higher. The prediction  accuracy of the Linear-SIF is significantly higher than that of the Fixed-SIF. When the loading angle is 15 • and 30 • , the Linear-SIF can achieve high precision prediction, and the maximum crack growth error is less than 1.8mm. When the loading angle is 45 • and 60 • , the maximum crack growth error is less than 4mm.
The prediction accuracy of the traditional FEM greatly depends on the crack growth increment a, and as the smaller value of a can achieve crack growth prediction performance. This is because the SIF under the current crack path can be calculated with a smaller a, but at the same time, the calculation efficiency problem will appear. Although the Linear-SIF has a higher prediction accuracy than the Fixed-SIF for the same a, its computational efficiency problem is more severe than the latter method because it requires two finite element simulations and linear interpolation of the SIF in each calculation cycle.
Compared with the traditional FEM, the fatigue crack growth prediction approach based on the improved GP mode can better fit the experimental data. Because the GP method can predict the change of the SIF caused by the crack growth at each step, and this change is also brought into the procedure for calculation. With the help of a high-performance computer, it not only improves the prediction accuracy of fatigue crack growth but also ensures a sufficiently fast calculation speed. Therefore, although the prediction results obtained by the GP method at 15 • and 30 • loading angles are slightly improved than the Linear-SIF, the calculation time is only 30% of the Linear-SIF with a = 0.6 mm (Calculation time does not include sample space establishment). This is because when the loading angle is small, the crack deflection angle will also be small, so the stress intensity factor at the crack tip does not change much. Therefore, when the loading angle is small, the main advantage of the approach in this paper is more reflected in the improvement of computational efficiency.
At 45 • and 60 • , the SIF at the tip of the crack changes more drastically due to the increase of the crack deflection angle. Therefore, the superiority of the calculation accuracy of the approach in this paper is reflected, and its calculation accuracy is higher than that of the Linear-SIF with a = 0.3 mm improves 39% and 61%. However, it is noted that the approach proposed in this paper predicts that the initial growth error of the crack at 60 • will reach a maximum of 0.7 mm, which also occurs because the SIF varies too much at small crack lengths, making the GP model unable to accurately predict such large changes. Therefore, the solution to this problem is to appropriately increase the number of sample points at the crack deflection point when the loading angle is large to allow the GP model to better predict this change.  The coordinate origin in the figure represents the tip point of the initial crack, while the x-axis and y-axis are the parallel and perpendicular directions of the initial crack, respectively. The final crack growth prediction path is almost the same as the prediction result using the GP method when the crack growth is a = 0.6 mm and a = 0.3 mm. According to the results, the difference between the GP method and the FEM method is small, and the error is within an acceptable range. This is because the total length of crack growth set for this fatigue crack growth experiment is small, and the experiment is stopped when the total length of crack growth under each load is only about 10 mm. However, in the CTS specimen, the small total length of crack growth leads to a small change in the deflection angle during crack growth prediction. Therefore, if the amount of crack growth is small, both of methods for predicting the crack growth path are acceptable. Fig. 18 shows the prediction of the crack growth path under various loading angles using the GP method. The error of the prediction results increases with the increase of the loading angle, but the accuracy of the path results predicted by the GP method is consistent with engineering demands.

VI. CONCLUSION
This paper proposes an approach for predicting fatigue crack growth based on the GP model with local sample densification and this new approach is suitable for crack growth prediction under a mixed-mode loading (Mode I+II). In this approach, the GP model is used instead of the traditional FEM to calculate the SIF at the crack tip, which realizes the nonlinear continuous change of this parameter in the fatigue crack growth prediction. Then, based on the assumptions of the crack growth path, an improved local sample densification approach is used to dynamically update the GP samples during the fatigue crack growth process, so that the number of samples can be reduced while ensuring the accuracy of prediction. Finally, the effectiveness of the approach proposed in this paper is verified according to the crack growth experimental of the lower bainite steel (SCM435) material, and the results show that: (1) By comparing the SIF-GP model using the ARD-Matérn (5) kernel function with the SIF calculated by the FEM, the maximum relative error does not exceed 1.5%, which proves that the GP model can solve the SIF for different crack paths with high accuracy. (2) By comparing the crack growth life predicted by the approach proposed in this paper and the traditional FEM, it is found that when the loading angle is small, the accuracy of the FEM is not much different from the approach based on GP model, but it can effectively improve the computational efficiency. When the loading angle is large, the approach based on GP model is not only more efficient but also more accurate than the FEM. (3) For the comparison between the two approaches for crack growth path prediction, it is found that when the total length of crack growth is small, the prediction accuracy between the approach based on GP model and the FEM is not much different, and the error will increase as the load loading angle increases. In the future, the following work will be performed: (1) it is necessary to compare the accuracy of the crack growth paths and to investigate the predictive ability of the method for fatigue crack growth under variable amplitude loads through relevant tests; (2) based on the current work, the proposed approach should be improved to make it suit to the prediction of fatigue crack growth in the mode-III.