Hybrid residual fatigue life prediction approach for gear based on Paris law and particle filter with prior crack growth information

2.1. Short overview of Paris law

In 1963, Paris, et al [19] proposed the Paris law based on fracture mechanics, which can reflect the failure mechanism of materials and is usually applied as a method to predict fatigue life or residual fatigue life. The stress intensity factor (SIF) plays a key role in the process of the fatigue crack propagation, and a number of tests are carried out to study the influence for crack propagation. It has been found that for the first loading mode the SIF increases along with the increase of the crack depth and for the second loading mode the SIF decreases along with the increase of the crack depth. The crack propagation rate $d\alpha /dN$ measured was similar with the SIF amplitude variation. Based on the above, the Paris law for the fatigue crack propagation under the constant amplitude loading was proposed:

where $\alpha$ represents the crack size, $N$ represents the number of stress cycles, $\mathrm{\Delta }K$ is the range of the stress intensity factor, $\mathrm{\Delta }\sigma$ is the stress range, $C$ and $m$ are the parameters of the Paris law which are usually determined by material constants.

2.2. Particle filter model

The particle filter employs the Monte Carlo simulation of a state dynamic model and Bayesian estimation for estimating the posterior probability density function (PDF) of the state. Therefore, the method is an effective tool for nonlinear and non-Gaussian systems. The state equation and measurement equation of particle filter are expressed as follows [20]:

where $x_k$ is the damage state to be estimated, $k$ is the time step index, ${\theta }_{k-1}$ is a vector of model parameters,$y_k$ is the measurement data, $n_k$ and ${\omega }_k$ are process and measurement noise, respectively. $f$ and $h$ represent the known process and observation functions, respectively.

The probability density function of the update unknown parameters can be obtained based on the following Bayes’ theorem:

where $\mathrm{\Theta }$ is a set of the unknown parameters, $z$ is a vector of monitoring data, $p\left({\mathrm{\Theta }}_{prior}\right)$ is the prior PDF of parameters, $p\left(\mathrm{\Theta }|z\right)$ is the posterior PDF of parameters based on the observations and $L\left(z|{\mathrm{\Theta }}_{prior}\right)$ is the likelihood based on the given parameters ${\mathrm{\Theta }}_{prior}$.

When the initial PDF of the parameters is given, in order to obtain the posterior probability density function $p\left({\theta }_k|z_{0:k}\right)$ of the unknown parameters, the state equation is defined as [21]:

where the notation 0:$k-1$ means a set of parameters from cycles 0 to $k-1$.

The new observations $z_k$ are collected at time $k$. According to the Bayesian rule, the posterior probability density for unknown parameters updates over time, so posterior probability distribution of the current parameters is obtained [21]:

where $p\left(z_k|z_{0:k-1}\right)=\int p\left({\theta }_k|z_{0:k-1}\right)p\left(z_k|{\theta }_k\right)d{\theta }_k$ is a constant.

Above is the prediction process and update process for parameters of Paris model. In this paper, we utilize the Bayesian filtering and Monte Carlo algorithm to obtain the parameters updating.

Assuming that update observations $z_{0:k}$ are known, the posterior probability density of parameters ${\theta }_{0:k}$ can be expressed as follows [21]:

where $\delta \left(\cdot \right)$ is the Dirac delta measure, ${\xi }_{0:k}$ is the priori vectors for parameters, ${\theta }_{0:k}$ is a set of vectors for parameters at cycles from 0 to k and $p\left({\xi }_{0:k}|z_{0:k}\right)$ is the priori probability density function. If we get the real posterior probability density function $p\left({\theta }_{0:k}|z_{0:k}\right)$, Eq. (6) can be calculated according to Eq. (7) [21]:

where ${\theta }_{0:k}^i,i=\mathrm{1,2}...,N_s$ are independent random samples sampling from $p\left({\theta }_{0:k}|z_{0:k}\right)$, $N_s$ is the number of samples sampling from $p\left({\theta }_{0:k}|z_{0:k}\right)$.

In the calculation process, $p\left({\theta }_{0:k}|z_{0:k}\right)$ are multi-type and non-standard. Therefore, it is difficult to obtain the result. In order to calculate the posterior probability density function of parameters for Paris model, the importance function sampling method is proposed in this paper. In this way, we can obtain the PDF distribution $\pi \left({\theta }_{0:k}|z_{0:k}\right)$ based on importance sampling for $p\left({\theta }_{0:k}|z_{0:k}\right)$, which has the same distribution with the $p\left({\theta }_{0:k}|z_{0:k}\right)$. Therefore, Eq. (7) can be transformed as follows [21]:

where $w_k^{*i}=\frac{p\left(z_{0:k}|{\theta }_{0:k}^i\right)p\left({\theta }_{0:k}^i\right)}{p\left(z_{0:k}\right)\pi \left({\theta }_{0:k}^i|z_{0:k}\right)}$ is the weights of the unknown parameters, ${\theta }_{0:k}^i$, $i=\mathrm{1,2},...,N_s$, which can be calculated by $\pi \left({\theta }_{0:k}^i|z_{0:k}\right)$sampling, and $p\left(z_{0:k}|{\theta }_{0:k}^i\right)$ is the probability density function of observations.

In order to obtain the weight, we assume that $p\left(z_{0:k}\right)=\int p\left(z_{0:k}|{\theta }_{0:k_{}}\right)p\left({\theta }_{0:k}\right)d{\theta }_{0:k}$, and the probability distribution $p\left({\theta }_{0:k}|z_{0:k}\right)$ can be expressed as follows:

where ${\tilde{w}}_k^{*i}=\frac{w_k^i}{\sum _{j=1}^{N_s}{w}_k^j}$ and $w_k^i=\frac{p\left(z_{0:k}|{\theta }_{0:k}^i\right)p\left({\theta }_{0:k}^i\right)}{\pi \left({\theta }_{0:k}^i|z_{0:k}\right)}=w_k^{*i}p\left(z_{0:k}\right)$.

3.1. Improved Paris law based on PF

In order to make full use of the advantages of data-driven method and model-based method, the integration of Paris law and PF model is proposed in this paper. The approach considers both the gear degradation process, and the prior information, which can describe the gear degradation process based on a simple model and improve the residual fatigue life prediction accuracy based on the prior information, which is the measured gear crack size as shown in Table 1. Based on the introduction of Paris law and PF model above, the state transition equation function for the Paris law is established, which is shown as follows:

The model parameters $m_k$ and $C_k$ as well as the degradation state ${\alpha }_k$ are estimated using the gear crack growth information $z_k$ under constant amplitude loading.

When the state transition equation and crack growth prior information are defined, the parameters $m_k$ and $C_k$ over time can be estimated based on the PF model. Moreover, we can estimate the parameters online when new gear crack growth information is added. In this paper, the parameters of the improved Paris law are obtained by maximum likelihood estimation, and the likelihood function is shown as follows:

where ${\zeta }_k^i=\sqrt{\mathrm{l}\mathrm{n}\left[1+{\left(\raisebox{1ex}{$\sigma $}\!\left/ \!\raisebox{-1ex}{${\alpha }_k^i\left(m_k^i,C_k^i\right)$}\right.\right)}^2\right]}$ and ${\lambda }_k^i=\mathrm{l}\mathrm{n}\left[{\alpha }_k^i\left(m_k^i,C_k^i\right)\right]-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\left({\zeta }_k^i\right)}^2$.

3.2. Implementation steps of hybrid prediction approach

In the prediction strategy process, the integration of the Paris model and particle filter is proposed to predict the residual fatigue life of the gear with prior crack growth information. The implementation steps of the proposed approach are shown as Fig. 1. The steps are described as follows:

Step 1. Decide the model-based prognostics method according to the object of study and raw measurement data. We choose the Paris law as the model-based method in this paper.

Step 2. Define the parameters which characterize the damage behavior. The model parameters which have an effect on the model behavior are often unknown and need to be identified in the prediction process.

Step 3. PF is employed to update the parameters as the data-driven model when new observations appear.

Step 4. If there is a new observation, go back to Step 3, or else go to Step 5.

Step 5. The residual fatigue life at the current time can be obtained based on the updated parameters of the Paris law.

The proposed hybrid approach considers the model-based prognostics and data-driven model at the same time, which can improve the prediction performance.

Fig. 1Implementation of hybrid prediction approach

4.1. Gear crack growth test

In order to verify the proposed method in this paper, a gear crack growth test is taken as an example to predict the residual fatigue life. As shown in Fig. 2, there is a through-the-thickness crack on the gear. The crack size is measured at every 50 cycles under the loading condition $\mathrm{\Delta }\sigma =$ 78 MPa, which is shown in Table 1, and the critical threshold of the crack size is 0.0463 m. The actual fatigue life is 2500 cycles. All the original data is taken from Ref [20]. In the current paper, we establish the model based on the crack size, not considering the relationship between the gear crack and condition monitoring data, such as vibration and temperature. Firstly, the true crack size data are generated according to Eq. (10). The measured crack size data are then generated by multiplying noise, which is lognormally distributed with standard deviation of 0.001/$a_k$ (m). Actually, it has been shown that the distribution of crack size follows a lognormal distribution [22].

Fig. 2Gear crack

4.2. Degradation process and parameters updating

As shown in Fig. 3, the measured crack size increases gradually over cycles, which means that the trend of the gear degradation. It is assumed that the standard deviation of measurement is known. Also, the initial distribution of the parameters and the likelihood function are, respectively, normal and lognormal distributions, which are as follows: ${\alpha }_0~N$(0.01, (5×10^–4)²), $m_0~N$(4, 0.2²), $\mathrm{l}\mathrm{o}\mathrm{g}\left(C_0\right)~N$(–22.33, 1.12²).

Fig. 3Monitoring crack size

According to the prior information of gear crack size data, we can predict the gear degradation process based on the proposed hybrid approach with the initial distribution of the parameters. The state transition for gear crack growth is given by Eq. (10).

The gear degradation in the future can be predicted, and the predicted crack size is shown in Table 2. Moreover, the gear crack growth is shown in Fig. 4, and we can tell that the gear is close to failure at 2500 cycles according to the failure threshold, which is fit with the actual residual life of the gear.

Fig. 4Crack growth prediction

While the gear degradation process is predicted based on the PF model, the parameters of Paris model are updated simultaneously. When the gear crack size in the future is predicted, we can update the model parameters based on the new obtained crack size data. Taking the parameters estimation at 1200 cycles and 1500 cycles as examples, initial value $m=$ 3.8, $C=$ 1.5×10^–10 and historical data are both generated based on PF to get the best estimation for the gear degradation state (${\alpha }_{1200}$, $m_{1200}$, $C_{1200}$) and (${\alpha }_{1500}$, $m_{1500}$, $C_{1500}$), which are shown in Fig. 5 and Fig. 6, and the model parameters are convergent to the optimum value quickly. According to the gear crack size and the model parameters, we can estimate the residual fatigue life at the current time.

Fig. 5Updated parameters at 1200 cycles

a)

b)

Fig. 6Updated parameters at 1500 cycles

a)

b)

4.3. Prediction results and discussion

The process of residual fatigue life prediction is shown in Fig. 7. Firstly, the prediction model needs to be defined based on the investigated object, and the model parameters are obtained according to the prior information, which has been finished in section 4.2. However, we can only get a lower accuracy result based on the model. Therefore, new monitoring information is added to the predicted process, and posterior parameters of the model can be obtained, which can improve the accuracy and update the parameters online. Finally, the residual fatigue life can be obtained according to the given failure threshold.

Fig. 7Illustration of RFL prediction

The model parameters are given at 1200 cycles and 1500 cycles as mentioned above. Consequently, the residual fatigue life distribution is shown in Fig. 8 and Fig. 9 respectively, which takes $N_s=$ 5000 particles in every simulation. According to the residual fatigue life distribution and the result of the residual fatigue life at the current time, we can obtain the probability density function, which is more suitable for describing the gear residual fatigue life.

Taking the calculation process above repeated, then the gear residual fatigue life can be obtained over time gradually. The gear residual fatigue life and 90 % confidence interval at a different time can be calculated based on the model when the parameters and failure threshold are given as shown in Table 3. The probability density function of the gear residual fatigue life is shown in Fig. 10, and the red curve and green curve represents the predicted mean residual fatigue life and actual residual fatigue life, respectively. The variances of the probability density function are smaller over cycles with more monitored and predicted crack size data is used during the prediction process, and it means that the prediction results are more and more accurate.

Fig. 8Residual fatigue life distribution at 1200 cycles

Fig. 9Residual fatigue life distribution at 1500 cycles

Fig. 10Residual fatigue life probability density function

The predicted residual fatigue life and actual fatigue life are compared as shown in Fig. 11. The predicted results are closer to the actual residual gear fatigue life with the increase of monitored data. The gear degradation process allows to the Paris law, and the model parameters are updated by the PF, which can make full use of both the prior information of gear crack, and the posterior information. However, we only study the residual fatigue life prediction under the constant amplitude loading, and the case under variable amplitude loading needs a further study. The results errors are caused by the monitoring information errors and the randomness of the proposed model, however, the model accuracy satisfies with the most working condition. In addition, when the gear is close to a failure, the error is much bigger, which is caused by the change of the mechanical properties.

Fig. 11Comparison of actual and predicted RFL