4. Application to crack growth prognostics

The code can be modified easily by the users for various applications. In this section, an example of crack growth is used to explain how the MATLAB code BM.m can be modified.

### 4.1 Model definition: crack growth

In addition to the trace of samples shown in Figure 6, it is possible to plot the histogram of model parameters using the same samples. The following commands

When the true degradation model is known, it is possible to compare the predicted degradation with the true one. The following MATLAB script plots the median and confidence intervals of the predicted degradation along with the true

The degradation curves up to 50 weeks are shown in Figure 8 for the battery example. The true degradation with b ¼ 0:012 is shown with the black curve. The red curves show 5, 50 (median) and 95 percentiles of the predicted degradation, which are caused by signiLevel = 5 (line 14). The plot also shows the threshold (green line) and measurement data (blue asterisk marks). Based on the true model, the end of life of the battery is tEOL ¼ 29:72 weeks, and thus, the true RUL should

be tRUL ¼ 20:72 weeks. The prediction shows that the median of RUL is

plot the histograms of all model parameters. figure; % histogram of parameters

for j=1:p;

end

figure;

Figure 8.

20

subplot(1,p,j);

hist(samplResul(j,:),30);

Fault Detection, Diagnosis and Prognosis

degradation and the threshold.

degraTrue=exp(-0.012.\*time);

plot(time,degraTrue,'k'); hold on; plot(time(1:k1),measuData,'\*b'); plot(time(k1:end),degraPI(:,1),'–r'); plot(time(k1:end),degraPI(:,2:3),':r'); plot([0 time(end)],[thres thres],'g');

degraPI=prctile(degraPredi',perceValue)';

xlabel('weeks'); ylabel('Relative C/1 capacity');

Comparison of the predicted degradation with the true degradation.

In fatigue crack growth, the failure criterion is given in terms of the crack size. Therefore, it would be appropriate to use the size of crack as a degradation feature. In this case, the degradation feature monotonically increases, while it was monotonically decreased for the battery example. Assuming that a through-the-thickness center crack exists in an infinite plate under mode-I loading condition, the rate of fatigue crack growth can be expressed using the Paris-Erdogan model as

$$\frac{\text{d}a}{\text{d}N} = \text{C}(\Delta K)^m \tag{16}$$

where a is the half crack size, N is the number of cycles, m and C are model parameters, <sup>Δ</sup><sup>K</sup> <sup>¼</sup> <sup>Δ</sup><sup>σ</sup> ffiffiffiffiffi <sup>π</sup><sup>a</sup> <sup>p</sup> is the range of stress intensity factor, and <sup>Δ</sup><sup>σ</sup> is the stress range. It is assumed that time is the number of fatigue loading cycles. For the consistent notation, the crack size and cycles are replaced with ~y ¼ a and t ¼ N in the following explanation. Since the degradation model requires the crack size as a function of time and model parameters, Eq. (16) can be integrated to obtain the following degradation model:

$$\tilde{\mathbf{y}}(t;m,\mathbf{C}) = \left[ t\mathbf{C} \left( \mathbf{1} - \frac{m}{2} \right) \left( \Delta \sigma \sqrt{\pi} \right)^{m} + a\_0^{\mathbf{1} - \frac{m}{2}} \right]^{\frac{2}{2-m}} \tag{17}$$

The system is under fatigue loading with the range of stress being Δσ ¼ 75MPa at each cycle. It is assumed that the health monitoring is performed every 50 cycles to measure the crack size yk until the current time tCUR ¼ 1, 200 cycles, and the initial size of the crack is a<sup>0</sup> ¼ 0:01m. Similar to the battery degradation example, the measurement data are simulated by adding random noise to the true crack size. First, the true crack size data are generated at every 50 cycles using Eq. (17) with <sup>m</sup>true <sup>¼</sup> <sup>3</sup>:8 and <sup>C</sup>true <sup>¼</sup> <sup>1</sup>:<sup>5</sup> � <sup>10</sup>�10. The measured crack size data are then generated by adding Gaussian noise ε � N 0; s <sup>2</sup> ð Þ, s <sup>¼</sup> <sup>0</sup>:0005m to the true crack sizes. The measured crack size data are used to identify three model parameters, θ ¼ f g m; ln ð Þ C ; s . In the Paris-Erdogan model, the y-intercept C is very small but changes its magnitude by several orders. Therefore, it would be better to identify logarithm of C. For RUL calculation, the critical crack size is determined as 0.043 m.

For the Bayesian method, it is necessary to define the prior distribution and the likelihood function. In the battery example, it was assumed that noise in data follows a normal distribution. However, when the distribution type of measurement noise is unknown, it is possible that the likelihood function might be different from the true noise distribution. The same is true for the prior/initial distribution. Therefore, it would be a good exercise to study the effect of different distribution types by changing the MATLAB codes. In this example, the lognormal distribution is employed for the likelihood function as:

$$f\left(y\_k|\theta\right) = \frac{1}{y\_k \zeta\_k \sqrt{2\pi}} \exp\left[-\frac{1}{2} \left(\frac{\ln y\_k - \eta\_k}{\zeta\_k}\right)^2\right], \quad k = 1, \ldots, N\_{\text{data}}\tag{18}$$

where ζ<sup>k</sup> ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln 1 þ s=~ykÞ 2 <sup>r</sup> <sup>h</sup> � i and <sup>η</sup><sup>k</sup> <sup>¼</sup> ln <sup>~</sup>ykÞ � <sup>0</sup>:<sup>5</sup> <sup>ζ</sup><sup>k</sup> ð Þ<sup>2</sup> � are the standard

deviation and mean of a lognormal distribution, respectively. In the above equation, ~ykð Þθ is the model prediction from Eq. (17) at time tk with given model parameters m and C.

Also, the prior distribution of each parameter is assumed as a normal distribution as f mð Þ� <sup>N</sup> <sup>4</sup>; <sup>0</sup>:2<sup>2</sup> � �, <sup>f</sup>ð Þ� ln<sup>C</sup> <sup>N</sup> �22:33; <sup>0</sup>:5<sup>2</sup> � �, and f sðÞ� <sup>N</sup> <sup>5</sup> � <sup>10</sup>�4; <sup>1</sup> � <sup>10</sup>�<sup>4</sup> � �<sup>2</sup> � �. Therefore, the joint prior distribution can be obtained from the independence

assumption as

$$f\_{\;\;\theta}(\theta) = f(m) \times f(\ln \mathcal{C}) \times f(\mathfrak{s}) \tag{19}$$

Initially, a fatigue crack grows slowly and then grows rapidly just before becoming unstable. The crack growth model in Eq. (17) is only valid when the crack growth is stable. Normally the threshold ythreshold is set before the crack becomes unstable. When the crack becomes unstable, Eq. (17) yield a complex number. Therefore, the last line of the model definition code identifies if a prediction yields a complex number and converts it to a real number by ignoring the

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB

In addition to the problem definition part, the posterior distribution part also needs to be modified because instead of a uniform distribution, a normal distribution is used for the prior distribution. Also, lognormal distributions are used instead of normal distributions for the likelihood function. The posterior distribution part

zeta=sqrt(log(1+(s./degraModel(k)).^2)); eta=log(degraModel(k))-0.5\*zeta.^2;

Similar to the battery example, the MATLAB code can be used to plot the trace

of MCMC sampling and the histogram of model parameters and RUL. Using a similar code provided in Section 3.4, the degradation trend could also be plotted. For example, Figure 5 shows the predicted degradation with the true degradation. Even if the median (1553 cycles) has a relatively large error with the true RUL

An important information that was not discussed before is the correlation between model parameters. It is well known that the Paris-Erdogan model parameters, m and C, are strongly correlated [20]. Therefore, it would be beneficial to plot the MCMC samples in the parameter space. The following MATLAB script plots the

[,post(i,j)]=BMappl(para,ParamName,time(1:k1),measuData,prioDisPar);

plot(samplResul(1,:),samplResul(2,:),'.','color',[0.5 0.5 0.5]); hold on

Figure 9 shows the MCMC samples in the parameter space along with the exact value of the parameters (star marker). The figure also shows the contour of the joint posterior PDF based on the grid method. It can be observed that the MCMC samples represent the joint posterior PDF well, and the joint PDF covers the true parameter values. However, due to a strong correlation between the two Paris-Erdogan model parameters, the joint PDF shows a narrow but long tail. Any combination of model parameters along the correlation line may yield a similar

(1709 cycles), the 95% confidence interval covers the true RUL.

[m, C]=meshgrid((3.4:0.009:4.3)',(-24:0.025:-21.5));

prior=prod(normpdf(param,prioDisPar(:,1),prioDisPar(:,2)));

likel=lognpdf(measuData(k),eta,zeta).\*likel;

complex part.

likel=1;

end

4.3 Results

in lines 71–74 needs to be modified as follows:

DOI: http://dx.doi.org/10.5772/intechopen.82781

for k=1:length(measuData)

MCMC samples of the parameters.

for i=1:101; for j=1:101

contour(m,C,real(post));

end; end

damage growth trend.

23

clear; clc; load('Crack at 1200.mat');

plot(3.8,log(1.5E-10),'pk','markersize',14);

para=[m(i,j); C(i,j); 0.0005];

poste=likel.\*prior;

The posterior distribution can be obtained by multiplying the prior distribution in Eq. (19) with the likelihood function in Eq. (18).

#### 4.2 Modifying the code

For the crack growth example, the code in Appendix needs to be changed as follows. First, the problem definition part in lines 2–15 is replaced with the following code:

%===== PROBLEM DEFINITION 1 (Required Parameters)================== WorkName='Crack'; %work results are saved by WorkName TimeUnit='cycles'; % time unit name time=(0:50:3600)'; % time including both at measurement and at prediction measuData=[1.03 1.00 0.96 1.14 1.13 1.10 1.15 1.15 1.19 1.19 1.14 ... 1.14 1.20 1.22 1.37 1.21 1.25 1.25 1.36 1.30 1.32 1.48 ... 1.52 1.47 1.59]'\*0.01; thres=0.043; % threshold - critical value ParamName=['m'; 'C'; 's']; % model parameters' name to be estimated prioDisPar=[4 0.2; -22.33 .5; 5E-4 1E-4]; % parameter prior distributions Ns=10000; % num. of samples for MCMC simulation burnIn=0.2; % burn-in fraction signiLevel=2.5; % significance level for C.I. and P.I. %===============================================================

A total of Ndata ¼ 25 measurement data are provided up to tCUR ¼ 1, 200 cycles, and the degradation is predicted until tend ¼ 3, 600 cycles. Since the prior distributions are assumed a normal distribution, the first column of prioDisPar is the mean, and the second column is the standard deviation.

Next, the model definition part in lines 65–67 is replaced with the following code:

```
%===== PROBLEM DEFINITION 2 (model equation)=======================
a0=0.01; dsig=75; coef=1-m/2;
degraModel=(t.*exp(C).*coef.*(dsig*sqrt(pi)).^m + a0.^coef).^(1./coef);
loca=imag(degraModel)�=0; degraModel(loca)=real(degraModel(loca));
%===============================================================
```
Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB DOI: http://dx.doi.org/10.5772/intechopen.82781

Initially, a fatigue crack grows slowly and then grows rapidly just before becoming unstable. The crack growth model in Eq. (17) is only valid when the crack growth is stable. Normally the threshold ythreshold is set before the crack becomes unstable. When the crack becomes unstable, Eq. (17) yield a complex number. Therefore, the last line of the model definition code identifies if a prediction yields a complex number and converts it to a real number by ignoring the complex part.

In addition to the problem definition part, the posterior distribution part also needs to be modified because instead of a uniform distribution, a normal distribution is used for the prior distribution. Also, lognormal distributions are used instead of normal distributions for the likelihood function. The posterior distribution part in lines 71–74 needs to be modified as follows:

```
prior=prod(normpdf(param,prioDisPar(:,1),prioDisPar(:,2)));
likel=1;
for k=1:length(measuData)
 zeta=sqrt(log(1+(s./degraModel(k)).^2)); eta=log(degraModel(k))-0.5*zeta.^2;
 likel=lognpdf(measuData(k),eta,zeta).*likel;
end
poste=likel.*prior;
```