- Email: [email protected]
Detecting cracks in aircraft engine fan blades using vibrothermography nondestructive evaluation
Reliability Engineering and System Safety 131 (2014) 229–235
Contents lists available at ScienceDirect
Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress
Detecting cracks in aircraft engine fan blades using vibrothermography nondestructive evaluation$ Chunwang Gao a, William Q. Meeker a,n, Donna Mayton b a b
Department of Statistics, Iowa State University, Ames, IA 50011, United States 21950 Nicole Street, Ramona, CA 92065, United States
art ic l e i nf o
a b s t r a c t
Available online 10 June 2014
Inspection is an important part of many maintenance processes, especially for safety-critical system components. This work was motivated by the need to develop more effective methods to detect cracks in rotating components of aircraft engines. This paper describes the analysis of data from vibrothermography inspections on aircraft engine turbine blades. Separate but similar analysis were done for two different purposes. In both analyses, we fit statistical models with random effects to describe the crackto-crack variability and the effect that the experimental variables have on the responses. In the first analysis, the purpose of the study was to find vibrothermography equipment settings that will provide good crack detection capability over the population of similar cracks in the particular kind of aircraft engine turbine blades that were inspected. Then, the fitted model was used to determine the test conditions where the probability of detection (POD) is expected to be high and probability of alarm is expected to be low. In our second analysis, crack size information was added and a similar model was fit. This model provides an estimate of POD as a function of crack size for specified test conditions. This function is needed as an input to models for planning in-service inspection intervals. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Probability of detection Random effects Safe life Thermal acoustics Thermosonics
1. Introduction 1.1. Vibrothermography in nondestructive evaluation Inspection is an important part of many maintenance processes for safety-critical system components, especially when using the damage tolerance approach to assure reliability, described for example, in Xiong and Shenoi [21]. Better inspection methods increase system reliability and are more economical because they allow longer periods of time between inspections. Because inspections are imperfect, it is critical to have an assessment of probability of detection (POD). See, for example, Badia et al. [2] and Cronvall et al. [4]. Simola and Pulkkinen [19] describe statistical methods for estimating POD from inspection data. Vibrothermography, also known as Thermal Acoustics (TA), thermosonics and Sonic Infrared (IR), is a relatively new technique used for nondestructive evaluation to detect fatigue cracks in metals and composites. An ultrasonic gun, which vibrates at a frequency range of 20 kHz–40 kHz, is used as the energy source.
☆ n
This paper is submitted to the RESS Special Issue on Accelerated Testing. Corresponding author. Tel.: þ 1 515 294 5336. E-mail address: [email protected] (W.Q. Meeker).
http://dx.doi.org/10.1016/j.ress.2014.05.009 0951-8320/& 2014 Elsevier Ltd. All rights reserved.
A short pulse (50–200 ms) of high power acoustical energy is applied to the test sample to make it vibrate. If a crack exists in the material, it is expected that the faces of the crack will rub against each other, resulting in a temperature increase near the crack. An infrared camera is used to record the temperature change in a sequence of images over time, which we call a movie. See, for example, Henneke and Jones [10], Shull [18], Hellier [9], Mix [15] and Favro et al. [6] for more details. 1.2. Purpose The purpose of the study described in this paper is to find vibrothermography equipment settings that will provide good crack detection capability over the population of cracks that could exist in aircraft engine turbine blades and to estimate the POD for specified test conditions as a function of crack size. We extend previous work of Mayton et al. [13] who modeled data from an experiment on turbine disks. This paper, however, is more general. We describe statistical methods that can be used for the same purpose in other inspection applications where there are equipment tuning parameters that must be set in order to have good inspection capability. In our study of turbine fan blades, there was an important amount of blade-to-blade variability that required us to us a mixed
230
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
effect regression model. We use the model to determine the test conditions where the probability of detection is expected to be high and probability of a test-stopping alarm is expected to be low. We fit a similar model that includes crack size as an explanatory variable and use this model to quantify Probability of Detection (POD) as a function of crack size. 1.3. Overview The remainder of this paper is organized as follows. Section 2 describes the experimental setup with data description. Section 3 builds a mixed effect model for the response variable. Probability of detection is studied in Section 4 and the alarm condition is studied in Section 5. Based on the result from Sections 4 and 5, a procedure to choose a good test condition is discussed in Section 6. Section 7 adds crack size information into the mixed effect model and the corresponding probability of detection methods is studied in Section 8. Also included in Section 8 are the confidence interval calculation and an evaluation of the coverage probability of the confidence intervals.
Table 1 Test factor levels. Level
Vibration amplitude (% of maximum)
Pulse length (ms)
Trigger force (pounds)
Low Medium High
10 35 60
83 150 217
10 35 60
Table 2 Sample test data of the crack in Blade 1, which include the test conditions and response results. Vibration amplitude
Pulse length
Trigger force
Maximum contrast
35 35 35
150 150 150
35 35 35
7 7 10
60 60 60
217 217 217
35 35 35
117 122 109
2. Experimental setup The energy for the vibrothermography system was generated by 2KW Branson model 921AES ultrasonic welder system. Data were acquired using an Indigo Merlin-Mid IR camera with a halfinch extender ring that was operated at 60 frames per second. A customized version of the EchoTherm s VT program, from Thermal Wave Imaging, Inc (TWI), was used to control both the IR camera and the ultrasonic welder. More details on the experimental set up can be found in Mayton et al. [13]. Fig. 1 shows a picture of a turbine blade used in this study. The test matrix consisted of experimental factors vibration amplitude, pulse length, and trigger force, each at three levels, as shown in Table 1. This test matrix was used to inspect 10 cracks in seven blades (two cracks were tested on Blade 4 and three cracks were tested on Blade 7). There were 32 tests (all 27 combinations of the 33 full factorial plus one additional test at 10-83-10, two additional tests at 35-150-35, and two additional tests at 60-217-60) for each crack. The response was maximum contrast, defined as the difference between peak value at the crack and background (3 3 pixel area 15 pixels to the left of the peak), were manually extracted. The responses (maximum contrasts) from the three replications in the same test were averaged before being used in model fitting. Table 2 gives sample data (two tests result) from the crack in Blade 1.
3. Analysis of signal response data 3.1. Transformation and fixed effects linear model Different model forms and transformation of variables were tried to find a model that adequately describes the relationship between the response and the test conditions. A linear model (i.e., a model for a possibly transformed response and explanatory
Table 3 Fitting result of model (4) for individual cracks. Crack b
β0;b
β1;b
β2;b
β3;b
1 2 3 4 5 6 7 8 9 10
6.580 20.505 38.022 12.915 27.055 6.967 11.945 37.345 4.654 179.016
7.651 4.909 5.006 3.180 7.914 0.326 3.160 0.495 1.119 1.227
2.068 8.423 11.639 4.118 15.356 5.131 3.948 10.730 3.567 59.186
5.428 0.287 4.157 1.793 1.762 2.438 0.155 6.463 0.153 23.576
variables that is linear in the parameters) with a log transformation on the response and the explanatory variables was chosen to describe the each individual crack. The linear model for the transformed maximum contrast response variable (denoted by M) can be expressed as log ðM þ1Þ ¼ β0;b þ β1;b log ðvÞ þ β2;b log ðpÞ þβ3;b log ðrÞ þ β4;b log ðvÞlog ðrÞ þ β5;b log ðvÞlog ðpÞ þ β6;b log ðpÞlog ðrÞ þ β7;b log ðvÞ2 þ β8;b log ðpÞ2 þ β9;b log ðrÞ2 þ ε
ð1Þ
where v represents Vibration Amplitude, p represents Pulse Length, and r represents Trigger Force. The effect of the þ1 in the transformation of the maximum response M is to weaken the effect of the log transformation, needed because there were some responses that were close to 0. The function lm in R [17] was used to fit the model using ordinary least squares. The estimates of the first four coefficients are listed in Table 3. The large crack to crack variability in these coefficients suggested a random effect model in which some of the parameters are random. 3.2. Full mixed effects model
Fig. 1. Aircraft engine turbine blade, approximately 5.5 inches long.
In order to describe a larger population of cracks using one model, we use a linear mixed effects model (described, for example, in Venables and Ripley [20] and Pinheiro and Bates [16]) in which the intercept and main effect (linear) model terms are modeled as random effects to describe the crack-to-crack variability in the response. We assume that the 10 cracks in the
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 4 Fixed effect estimates of the full mixed-effect model.
231
Table 6 Fixed effect estimates for the reduced mixed-effect model.
βs
Value
Std. error
t-ratio
p-value
βs
Value
Std. error
t-ratio
p-value
β0 β1 β2 β3 β4 β5 β6 β7 β8 β9
7.843 3.580 2.893 1.137 0.671 0.092 0.299 0.214 0.265 0.129
11.620 1.272 4.672 1.427 0.166 0.086 0.173 0.152 0.472 0.163
0.675 2.814 0.619 0.797 4.034 1.061 1.729 1.409 0.561 0.791
0.500 0.005 0.536 0.426 0.000 0.290 0.085 0.160 0.575 0.429
β0 β1 β5 β9
2.974 1.248 0.521 0.139
0.327 0.442 0.076 0.026
9.085 2.822 6.858 5.356
0.000 0.005 0.000 0.000
Table 7 Covariance matrix for the random effects for the reduced mixed effect model.
Table 5 Covariance matrix of the random effect model. Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
19.073 1.004 3.688 2.120
1.004 0.117 0.215 0.126
3.688 0.215 0.721 0.421
2.120 0.126 0.421 0.281
study represent a random sample from that population of cracks. The final model for the maximum contrast is
Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
18.079 0.935 3.482 2.034
0.935 0.109 0.199 0.119
3.482 0.199 0.678 0.403
2.034 0.119 0.403 0.273
The model assumptions can be stated succinctly as ε Normal ð0; σ ε Þ and 0 1 00 1 1 β0;b β0 Bβ C BB β C C B 1;b C BB 1 C C B C B C; Σ βb C: B β2;b C Multinormal B @ A @ A 0 @ A β3;b 0
log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ 4. Probability of detection for a population of cracks
β4 log ðvÞlog ðrÞ þ β5 log ðvÞlog ðpÞ þ β6 log ðpÞlog ðrÞ þ β7 log ðvÞ2 þ β8 log ðpÞ2 þβ9 log ðrÞ2 þ ε:
ð2Þ
The random error ε is assumed to have a normal distribution with mean 0 and variance σ ε . For each crack, the coefficients ðβ0;b ; β1;b ; β2;b ; β3;b Þ are assumed to have a multivariate normal distribution independent of ε, with mean μβb ¼ ðβ0 ; β1 ; β2 ; β3 Þ and covariance matrix Σ βb , where μβb is a vector and Σ βb is a 4 4 matrix. These model assumptions can be stated succinctly as ε Normalð0; σ ε Þ and 0 1 00 1 1 β0;b β0 Bβ C BB C C B 1;b C BB β C C B C Multinormal BB 1 C; Σ βb C: B β2;b C BB β 2 C C @ A @@ A A β3;b β3 The parameter estimates, which are calculated by the function lme (in package nlme) in R [17], are given in Table 4 and the matrix elements estimates are given in Table 5.
3.3. Reduced mixed effect model We reduced the full model in (2) based on Akaike's information criterion (AIC) [1] and the t ratios of the regression coefficients. The final model that we used is log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β5 log ðvÞlog ðpÞ þ β9 log ðrÞ2 þ ε:
ð3Þ
We still assumed that ðβ0;b ; β1;b ; β2;b ; β3;b Þ has a multivariate normal distribution independent of ε, but with fixed 0 mean for β2;b and β3;b , i.e., μβb ¼ ðβ0 ; β1 ; 0; 0Þ. The reduced model sets β5, β7 and β8 to be 0. The results for fitting the reduced model are listed in Table 6. The covariance matrix results are given in Table 7.
4.1. Distribution of response variable Under the reduced mixed effect model (3), for given fixed values of v, p and r, the response variable m þ1 has a normal distribution with mean μðv; p; rÞ ¼ β0;b þ β1;b log ðvÞ þ β2;b log ðpÞ þβ3;b log ðrÞ þ β5 log ðvÞlog ðpÞ þ β9 log ðrÞ2
ð4Þ
and variance σ 2 ðv; p; rÞ ¼ X Tl Σ βb X l þ σ 2ϵ :
ð5Þ
Here X l ¼ ð1; log ðvÞ; log ðpÞ; log ðrÞÞT where superscript T indicates vector transpose. Eq. (5) will provide basis for estimating probability of detection (POD) for the population of potential cracks as a function of the experimental factors. 4.2. Calculation of probability of detection Probability of detection (POD) is described in detail in MILHDBK-1823A [14] and Georgiou [8]. In our application, POD is the probability that M from a thermal image sequence exceeds a threshold mT for a crack taken at random from the population of cracks. This probability can be expressed as log ðmT þ1Þ μðv; p; rÞ ð6Þ PODðmT ; v; p; rÞ ¼ PrðM 4mT Þ ¼ 1 Φ σ where ΦðxÞ is the standard normal cumulative distribution function distribution. We can estimate the POD by evaluating Eq. (6) using the estimates of the model parameters in Table 4. The threshold was chosen as mT ¼10. This value was chosen to balance sensitivity with the probability of generating false positives. Fig. 2 displays the output of our mixed-effects model for the POD. The darker the color, the lower the POD value.
232
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Fig. 2. POD plot for pulse length¼ 217. The lighter the color, the higher of the POD.
5. Analysis of the alarm data
Fig. 3. POA plot for pulse length¼ 217. The lighter the color, the larger of the POA.
Some tests ended with an alarm, usually caused by a lock-up of the welder/sample system when the horn can no longer vibrate properly. This is caused by protective feedback loop that terminates power to the welder to prevent damage to the piezoelectric element. Alarm incidents are undesirable because no useful information is obtained from the test. We would like to find test conditions that will also result in a small probability of an alarm. When there was an alarm, the time to alarm was recorded. When there was no alarm, the alarm time was taken to be right censored at the length of the pulse because the unknown alarm time is greater than the pulse length. We fit a log normal regression model to the censored time-to-alarm data to assess the probability of alarm as a function of the test conditions. Let TA denote the time of an alarm. A probability plot of the censored data suggested that TA can be described by a lognormal distribution so that log ðpÞ μa POA ¼ PrðT A rpÞ ¼ Φnor : ð7Þ σa Here μa ¼ βa0 þβa1 v þ βa2 r. Higher order terms were not statistically significant. The restricted maximum likelihood (REML) estimates a a a of the parameters are β^ 0 ¼ 6:027, β^ 1 ¼ 0:0188, β^ 2 ¼ 0:00531, a σ^ ¼ 0:1625. Fig. 3 depicts estimates of POA for pulse length equal to 217 ms.
Fig. 4. Overlay threshold plot for pulse length ¼217. The South-east corner of the plot satisfied both criteria (i.e., POD 4 0:90 and POA o 0:05).
6. Choosing good test conditions
7. Modeling response as a function of effective crack size
A choice of good test conditions needs to fulfill two requirements: a large response for the crack and a small chance of triggering an alarm that would terminate the test and yield no information. Fig. 4 combines the information in the POD plot in Fig. 2 and the POA plot in Fig. 3, showing the test conditions where the probability of detection is expected to be high and probability of alarm is expected to be low. Such overlay threshold plots require specification of acceptable POD and PFA thresholds and we chose those to be POD 4 0:90 and POA o0:05. The south-east region of Fig. 4 meets both conditions, suggesting the test conditions to be used in future tests.
In this section, we will build a statistical model to describe signal amplitude as a function of crack length. This model can then be used to estimate POD as a function of crack length. 7.1. Crack sizes Table 8 lists the crack size information that was obtained by using a method known as acetate replication, a nondestructive technique. Apparent length is the portion of the total wrap-around crack length visible from the same side of the airfoil as the IR camera. Edge length is the portion of the crack length that wraps
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 8 Size information for the 10 cracks.
233
Table 10 Covariance matrix of random effect of full model with size.
Crack
Apparent
Opposite apparent
Edge
Total
Effect
Σ βb
β0
β1
β2
β3
1 2 3 4A 4B 5 6 7A 7B(& C) 7D
0.04 0.1058 0.2903 0.0793 0.2907 0.0935 0 0.2781 0.0962 0.1039
0 0.0304 0.228 0.1875 0 0.1129 0.0718 0.2449 0.1321 0.0503
0.0039 0.021 0.0153 0.0313 0 0.0163 0.0066 0.0156 0.01 0.0126
0.0439 0.1572 0.5336 0.2981 0.2907 0.2227 0.0784 0.5386 0.2383 0.1668
0.0046 0.0226 0.0204 0.0322 0.0054 0.0177 0.0065 0.0205 0.0116 0.0143
β0 β1 β2 β3
21.158 0.929 3.797 2.336
0.929 0.114 0.206 0.128
3.797 0.206 0.705 0.434
2.336 0.128 0.434 0.270
Table 9 Fixed effect estimates of the mixed-effect model with size. βs
Value
Std. error
t-ratio
p-value
β0 β1 β2 β3 β4 β5 β6 β7 β8 β9 β10
3.095 3.517 2.874 1.185 0.662 0.087 0.300 0.213 0.260 0.124 1.122
11.765 1.283 4.713 1.440 0.168 0.087 0.175 0.153 0.476 0.164 0.187
0.263 2.741 0.610 0.823 3.944 0.996 1.716 1.391 0.546 0.755 6.006
0.793 0.007 0.543 0.412 0.000 0.320 0.088 0.165 0.586 0.451 0.000
Table 11 Fixed effect estimates of the reduced mixed-effect model with m as a function of crack size. βs
Value
Std. error
t-ratio
p-value
β0 β1 β5 β6 β9 β10
3.010 1.774 0.637 0.164 0.253 1.127
1.145 0.564 0.108 0.107 0.079 0.186
2.628 3.144 5.874 1.531 3.192 6.048
0.009 0.002 0.000 0.127 0.002 0.000
β3;b (i.e., μβb ¼ ðβ0 ; β1 ; 0; 0Þ). The results for fitting the reduced model are given in Table 11 and the covariance matrix results are given in Table 12.
8. POD as a function of effective crack size 8.1. POD versus crack size
around the edge of the airfoil but is not visible from either side of airfoil. Opposite apparent length is the portion of the crack that would be visible from the opposite side of the airfoil. Especially because there were zeros for some of these reported sizes, we use a convex combination of these recorded sizes to define an effective crack length. We then treat this variable as a fixed effect in the regression model. 7.2. Full mixed effect model with crack size information The full mixed effects model that includes crack length, shown in (8), is similar to model (2). log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β4 log ðvÞlog ðpÞ þ β5 log ðvÞlog ðrÞ þ β6 log ðpÞlog ðrÞ þ β7 log ðvÞ2 þ β8 log ðpÞ2 þβ9 log ðrÞ2 þ β10 log ðLN Þ þ ε:
ð8Þ
Here LN ¼ αLa þ ð1 αÞLe is the effective crack length and α¼0.0185 was chosen to maximize the likelihood. That is, a sequence of models was fit using different values of α between 0 and 1. The value of α that gave the largest value of the likelihood was α¼ 0.0185. The fixed effect parameter estimates of model (8) are given in Table 9. The corresponding estimates of the variances and covariance of the random effect terms are given in Table 10.
It is necessary to estimate POD as a function of crack size for purposes such as planning periodic inspections. In this section, we are going to use model (9) to recalculate the POD function. In contrast to Eq. (6), POD will now also be a function of effective crack size and can be expressed as log ðmT þ 1Þ μðv; p; r; lÞ PODðm; v; p; r; lÞ ¼ Prðm 4 mT Þ ¼ 1 Φ : σ ð10Þ An estimate of POD as a function of effective crack size can be computed by substituting the model estimates from Tables 11 and 12 into model (9). Fig. 5 shows the POD estimate as a function of effective crack length for different test conditions (Vibration Amplitude, Pulse Length and Trigger Force, respectively). These particular conditions were chosen to provide examples of different POD values. The inside tic marks at the bottom of the plot show the measured effective sizes of the cracks in the data set. The POD estimates to the right of effective size 0.03 involve extrapolation. Both of the full model estimates and reduced model estimates are shown in the plot indicating little difference between the two models when estimating POD. In the following sections wee use the reduced model for its simplicity. 8.2. Bootstrap confidence interval of POD estimation
7.3. Reduced mixed effect model with crack size information Similar to what we did in Section 3.3, the full model can be reduced to a smaller model based on AIC. The model suggested by this criterion is log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β5 log ðvÞlog ðrÞ þ β6 log ðpÞlog ðrÞ þβ9 log ðrÞ2 þβ10 log ðLN Þ þ ε:
ð9Þ
We again assume that ðβ0;b ; β1;b ; β2;b ; β3;b Þ has a multivariate normal distribution independent of ε, but with fixed 0 mean for β2;b and
The complicated model fitting procedure and calculation of POD make the direct computation of a confidence interval of POD difficult. Bootstrap methods provide a convenient alternative for such complicated procedures as described in Efron and Tibshirani [5]. The bootstrap procedures that we used are fully parametric. That is, we simulate bootstrap samples from our assumed parametric model and the REML estimates [3]. Fig. 6 shows the result for both the percentile and the BCA (bias-corrected accelerated) pointwise nonparametric bootstrap confidence intervals for POD
234
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 12 Estimated covariance matrix of random effect for the reduced model with m as a function of size. Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
20.302 0.884 3.633 2.257
0.884 0.109 0.197 0.123
3.633 0.197 0.674 0.419
2.257 0.123 0.419 0.263
1.00
POD vs. Effective Crack Length Contrast Threshold = 10
Coverage Probability
0.95
● ●
● ●
0.90
●
●
●
●
●
●
●
●
0.85
● ● ●
0.80
●
●
●
●
BCA 95% Confidence Interval, 10 Samples Percentile 95% Confidence Interval, 10 Samples BCA 95% Confidence Interval, 20 Samples Percentile 95% Confidence Interval, 20 Samples BCA 95% Confidence Interval, 40 Samples Percentile 95% Confidence Interval, 40 Samples
0.75
0.8
0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040
Partial Model Full Model 0.6
● ●
1.0 (60,217,60)
●
Crack Length Fig. 7. Coverage probability of POD Confidence Intervals calculated by the BCA and the percentile bootstrap method for 10, 20 and 40 crack samples. The solid line is by the BCA method and the dotted line is by the percentile method. For the range of the crack size of our sample, the BCA method has better coverage than the percentile method.
POD
(35,150,35)
0.4
versus effective crack length in test condition (35, 150, 35). 5000 bootstrap samples were used in the calculation.
0.2
8.3. Coverage Probabilities of the Percentile and the BCA Confidence Interval
(10,83,10) 0.0 0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
Effective Crack Length
Fig. 5. POD as a function of effective crack length. The solid line is for the reduced model and the dotted line is for the full model. The three test condition parameters (Vibration Amplitude, Pulse Length and Trigger Force) are indicated in the parentheses near the lines.
1.0
Although bootstrap procedures are expected to provide actual coverage probabilities that are closer to nominal than simpler methods such as the Wald's approximate method, the actual coverage probability, in general, will only be approximate. We used simulations to compare the coverage probabilities of the two bootstrap confidence interval calculation methods: the sample percentile method and the more refined the BCA method using 2000 new data sets were simulated from the fitting result of the reduced model and the percentile and the BCA bootstrap methods were used to calculate confidence intervals. The observed coverage probability is defined as Coverage ¼
0.8
0.6 POD
BCA Confidence Interval Percentile Confidence Interval 0.4
0.2
0.0 0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
Effective Crack Length
Fig. 6. 95% confidence intervals based on the percentile and the BCA bootstrap methods for test condition (Vibration Amplitude ¼35, Pulse Length¼ 150 and Trigger Force¼35). The solid line in the middle is the POD estimate from the mixed effect model. The dotted lines are the pointwise confidence intervals, where the large dotted line is for the BCA method and the small dotted line is for the percentile method.
NI N
ð11Þ
where NI denotes the number of confidence intervals that include the true value of POD (which is calculated in the model from which the data sets are simulated) and N is the total number simulated data sets. The simulation results of the coverage probability are shown in Fig. 7, for data sets with 10, 20 and 40 experimental units with three different symbols. The solid lines are for the BCA method and the dotted lines are for the percentile method. The BCA shows closer to nominal converge probabilities than the percentile method, as expected. The coverage probability of both the BCA method and the percentile method become closer to nominal when the sample size increases. The BCA method with 40 experimental units is the closest to nominal among all six evaluations. When sample size is small, (i.e., 10), the coverage probabilities are far away from nominal for larger cracks. As the sample size increases, the difference in coverage between small cracks and large cracks becomes smaller.
9. Conclusion and future work This paper has shown how to find good vibrothermography inspection test conditions for detecting cracks in fan blades and
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
how to estimate POD as a function of the effective size of a crack. The same general approach could be use for other kinds of inspections where tuning parameters need to be set. There are a number of areas for possible future research in this area.
It would be useful to develop methods for quantifying and displaying uncertainty in the overlay threshold plots like Fig. 2.
In our work, the standard maximum contrast response was
used. It would be useful to explore the use of other response variables in an attempt to use more completely the information in the sequence of image produced by a vibrothermography test system. For example, Li et al. [11] used a signal-to-noise ratio criterion after a matched filter as a response from a different vibrothermography inspection experiment. Gao and Meeker [7] use a method based on principal components analysis of a set of coefficients from a series of fitted regression models to extract a scalar measure of increase in temperature from the sequence-or-image data. With the use of diffuse prior information, Bayesian methods will provide results that are similar to the REML method used in our study. The Bayesian method provides an attractive method to construct confidence intervals for complicated functions of the parameters (such as POD) and results in other areas of application have shown that intervals computed in this manner have good frequentist coverage properties. Li et al. [12] use Bayesian methods with a mixed effect model to estimate POD based on ultrasonic inspection of synthetic hard alpha inclusions in a disk forging.
Acknowledgments This material used in the paper is based upon work supported by the Air Force Research Laboratory under Contract # FA8650-04C-5228 at Iowa State University's Center for NDE.
235
References [1] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control 1974;19(6):716–23. [2] Badia FG, Berrade MD, Campos CA. Optimal inspection and preventive maintenance of units with revealed and unrevealed failures. Reliab Eng Syst Saf 2002;78:157–63. [3] Bartlett MS. Properties of sufficiency and statistical tests. Proc R Soc Lond Ser A Math Phys Sci 1937;160(901):268–82. [4] Cronvall O, Simola K, Männistö I, Gunnars J, Alverlind L, Dillström P, et al. A study on the effect of flaw detection probability assumptions on risk reduction achieved by non-destructive inspection. Reliab Eng Syst Saf 2012;105:90–6. [5] Efron B, Tibshirani RJ. An introduction to the bootstrap. Chapman & Hall/CRC; 1994. [6] Favro LD, Han X, Ouyang Z, Sun G, Sui H, Thomas RL. Infrared imaging of defects heated by a sonic pulse. Rev Sci Instr 2000;71(6):2418–21. [7] Gao C, Meeker WQ. A statistical method for crack detection using vibrothermography inspection data. Qual Technol Quant Manag 2012;9:59–77. [8] Georgiou GA. Probability of detection (pod) curves derivation, applications and limitations. Technical report, Jacobi Consulting Ltd, 2006. Downloaded from 〈http://www.hse.gov.uk/research/rrpdf/rr454.pdf〉 on May 6; 2011. [9] Hellier C. Handbook of nondestructive evaluation. New York: McGraw-Hill Professional; 2001. [10] Henneke EG, Jones TS. Detection of damage in composite materials by vibrothermography. ASTM Spec Tech Publ 1979;1(696):83–95. [11] Li M, Holland SD, Meeker WQ. Statistical methods for automatic crack detection based on vibrothermography sequence-of-images data (with discussion). Appl Stoch Models Bus Ind 2010;26:481–95. [12] Li M, Meeker WQ, Thompson RB. Physical model-assisted probability of detection of flaws in titanium forgings using ultrasonic nondestructive evaluation. Technometrics 2014;56:78–91. [13] Mayton D, Spencer F, Alvarez C. Characterizing the effects of sonic IR variables on turbine disk inspection using a design of experiments approach. Rev Progr Quant Nondestruct Eval 2005;24:642–9. [14] MIL-HDBK-1823A. Nondestructive evaluation system reliability assessment, Standardization Order Desk, Building 4D, 700 Roberts Avenue, Philadelphia, PA 19111-5094, 2009. Downloaded from 〈http://www.statisticalengineering. com/mh1823/MIL-HDBK-1823A%282009%29.pdf〉 on May 6; 2011. [15] Mix PE. Introduction to nondestructive testing: a training guide. New York: Wiley-Interscience; 2005. [16] Pinheiro JC, Bates DM. Mixed-effects models in S and S-Plus. New York: Springer; 2000. [17] R Development Core Team, R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2008. [18] Shull PJ. Nondestructive evaluation: theory, techniques, and applications. New York: Marcel Dekker, Inc.; 2002. [19] Simola K, Pulkkinen U. Models for non-destructive inspection data. Reliab Eng Syst Saf 1998;60:1–12. [20] Venables WN, Ripley BD. Modern applied statistics with S. Fourth edition. New York: Springer; 2002. [21] Xiong JJ, Shenoi RA. A durability model incorporating safe life methodology and damage tolerance approach to assess first inspection and maintenance period for structures. Reliab Eng Syst Saf 2009;94:1251–8.
Contents lists available at ScienceDirect
Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress
Detecting cracks in aircraft engine fan blades using vibrothermography nondestructive evaluation$ Chunwang Gao a, William Q. Meeker a,n, Donna Mayton b a b
Department of Statistics, Iowa State University, Ames, IA 50011, United States 21950 Nicole Street, Ramona, CA 92065, United States
art ic l e i nf o
a b s t r a c t
Available online 10 June 2014
Inspection is an important part of many maintenance processes, especially for safety-critical system components. This work was motivated by the need to develop more effective methods to detect cracks in rotating components of aircraft engines. This paper describes the analysis of data from vibrothermography inspections on aircraft engine turbine blades. Separate but similar analysis were done for two different purposes. In both analyses, we fit statistical models with random effects to describe the crackto-crack variability and the effect that the experimental variables have on the responses. In the first analysis, the purpose of the study was to find vibrothermography equipment settings that will provide good crack detection capability over the population of similar cracks in the particular kind of aircraft engine turbine blades that were inspected. Then, the fitted model was used to determine the test conditions where the probability of detection (POD) is expected to be high and probability of alarm is expected to be low. In our second analysis, crack size information was added and a similar model was fit. This model provides an estimate of POD as a function of crack size for specified test conditions. This function is needed as an input to models for planning in-service inspection intervals. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Probability of detection Random effects Safe life Thermal acoustics Thermosonics
1. Introduction 1.1. Vibrothermography in nondestructive evaluation Inspection is an important part of many maintenance processes for safety-critical system components, especially when using the damage tolerance approach to assure reliability, described for example, in Xiong and Shenoi [21]. Better inspection methods increase system reliability and are more economical because they allow longer periods of time between inspections. Because inspections are imperfect, it is critical to have an assessment of probability of detection (POD). See, for example, Badia et al. [2] and Cronvall et al. [4]. Simola and Pulkkinen [19] describe statistical methods for estimating POD from inspection data. Vibrothermography, also known as Thermal Acoustics (TA), thermosonics and Sonic Infrared (IR), is a relatively new technique used for nondestructive evaluation to detect fatigue cracks in metals and composites. An ultrasonic gun, which vibrates at a frequency range of 20 kHz–40 kHz, is used as the energy source.
☆ n
This paper is submitted to the RESS Special Issue on Accelerated Testing. Corresponding author. Tel.: þ 1 515 294 5336. E-mail address: [email protected] (W.Q. Meeker).
http://dx.doi.org/10.1016/j.ress.2014.05.009 0951-8320/& 2014 Elsevier Ltd. All rights reserved.
A short pulse (50–200 ms) of high power acoustical energy is applied to the test sample to make it vibrate. If a crack exists in the material, it is expected that the faces of the crack will rub against each other, resulting in a temperature increase near the crack. An infrared camera is used to record the temperature change in a sequence of images over time, which we call a movie. See, for example, Henneke and Jones [10], Shull [18], Hellier [9], Mix [15] and Favro et al. [6] for more details. 1.2. Purpose The purpose of the study described in this paper is to find vibrothermography equipment settings that will provide good crack detection capability over the population of cracks that could exist in aircraft engine turbine blades and to estimate the POD for specified test conditions as a function of crack size. We extend previous work of Mayton et al. [13] who modeled data from an experiment on turbine disks. This paper, however, is more general. We describe statistical methods that can be used for the same purpose in other inspection applications where there are equipment tuning parameters that must be set in order to have good inspection capability. In our study of turbine fan blades, there was an important amount of blade-to-blade variability that required us to us a mixed
230
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
effect regression model. We use the model to determine the test conditions where the probability of detection is expected to be high and probability of a test-stopping alarm is expected to be low. We fit a similar model that includes crack size as an explanatory variable and use this model to quantify Probability of Detection (POD) as a function of crack size. 1.3. Overview The remainder of this paper is organized as follows. Section 2 describes the experimental setup with data description. Section 3 builds a mixed effect model for the response variable. Probability of detection is studied in Section 4 and the alarm condition is studied in Section 5. Based on the result from Sections 4 and 5, a procedure to choose a good test condition is discussed in Section 6. Section 7 adds crack size information into the mixed effect model and the corresponding probability of detection methods is studied in Section 8. Also included in Section 8 are the confidence interval calculation and an evaluation of the coverage probability of the confidence intervals.
Table 1 Test factor levels. Level
Vibration amplitude (% of maximum)
Pulse length (ms)
Trigger force (pounds)
Low Medium High
10 35 60
83 150 217
10 35 60
Table 2 Sample test data of the crack in Blade 1, which include the test conditions and response results. Vibration amplitude
Pulse length
Trigger force
Maximum contrast
35 35 35
150 150 150
35 35 35
7 7 10
60 60 60
217 217 217
35 35 35
117 122 109
2. Experimental setup The energy for the vibrothermography system was generated by 2KW Branson model 921AES ultrasonic welder system. Data were acquired using an Indigo Merlin-Mid IR camera with a halfinch extender ring that was operated at 60 frames per second. A customized version of the EchoTherm s VT program, from Thermal Wave Imaging, Inc (TWI), was used to control both the IR camera and the ultrasonic welder. More details on the experimental set up can be found in Mayton et al. [13]. Fig. 1 shows a picture of a turbine blade used in this study. The test matrix consisted of experimental factors vibration amplitude, pulse length, and trigger force, each at three levels, as shown in Table 1. This test matrix was used to inspect 10 cracks in seven blades (two cracks were tested on Blade 4 and three cracks were tested on Blade 7). There were 32 tests (all 27 combinations of the 33 full factorial plus one additional test at 10-83-10, two additional tests at 35-150-35, and two additional tests at 60-217-60) for each crack. The response was maximum contrast, defined as the difference between peak value at the crack and background (3 3 pixel area 15 pixels to the left of the peak), were manually extracted. The responses (maximum contrasts) from the three replications in the same test were averaged before being used in model fitting. Table 2 gives sample data (two tests result) from the crack in Blade 1.
3. Analysis of signal response data 3.1. Transformation and fixed effects linear model Different model forms and transformation of variables were tried to find a model that adequately describes the relationship between the response and the test conditions. A linear model (i.e., a model for a possibly transformed response and explanatory
Table 3 Fitting result of model (4) for individual cracks. Crack b
β0;b
β1;b
β2;b
β3;b
1 2 3 4 5 6 7 8 9 10
6.580 20.505 38.022 12.915 27.055 6.967 11.945 37.345 4.654 179.016
7.651 4.909 5.006 3.180 7.914 0.326 3.160 0.495 1.119 1.227
2.068 8.423 11.639 4.118 15.356 5.131 3.948 10.730 3.567 59.186
5.428 0.287 4.157 1.793 1.762 2.438 0.155 6.463 0.153 23.576
variables that is linear in the parameters) with a log transformation on the response and the explanatory variables was chosen to describe the each individual crack. The linear model for the transformed maximum contrast response variable (denoted by M) can be expressed as log ðM þ1Þ ¼ β0;b þ β1;b log ðvÞ þ β2;b log ðpÞ þβ3;b log ðrÞ þ β4;b log ðvÞlog ðrÞ þ β5;b log ðvÞlog ðpÞ þ β6;b log ðpÞlog ðrÞ þ β7;b log ðvÞ2 þ β8;b log ðpÞ2 þ β9;b log ðrÞ2 þ ε
ð1Þ
where v represents Vibration Amplitude, p represents Pulse Length, and r represents Trigger Force. The effect of the þ1 in the transformation of the maximum response M is to weaken the effect of the log transformation, needed because there were some responses that were close to 0. The function lm in R [17] was used to fit the model using ordinary least squares. The estimates of the first four coefficients are listed in Table 3. The large crack to crack variability in these coefficients suggested a random effect model in which some of the parameters are random. 3.2. Full mixed effects model
Fig. 1. Aircraft engine turbine blade, approximately 5.5 inches long.
In order to describe a larger population of cracks using one model, we use a linear mixed effects model (described, for example, in Venables and Ripley [20] and Pinheiro and Bates [16]) in which the intercept and main effect (linear) model terms are modeled as random effects to describe the crack-to-crack variability in the response. We assume that the 10 cracks in the
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 4 Fixed effect estimates of the full mixed-effect model.
231
Table 6 Fixed effect estimates for the reduced mixed-effect model.
βs
Value
Std. error
t-ratio
p-value
βs
Value
Std. error
t-ratio
p-value
β0 β1 β2 β3 β4 β5 β6 β7 β8 β9
7.843 3.580 2.893 1.137 0.671 0.092 0.299 0.214 0.265 0.129
11.620 1.272 4.672 1.427 0.166 0.086 0.173 0.152 0.472 0.163
0.675 2.814 0.619 0.797 4.034 1.061 1.729 1.409 0.561 0.791
0.500 0.005 0.536 0.426 0.000 0.290 0.085 0.160 0.575 0.429
β0 β1 β5 β9
2.974 1.248 0.521 0.139
0.327 0.442 0.076 0.026
9.085 2.822 6.858 5.356
0.000 0.005 0.000 0.000
Table 7 Covariance matrix for the random effects for the reduced mixed effect model.
Table 5 Covariance matrix of the random effect model. Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
19.073 1.004 3.688 2.120
1.004 0.117 0.215 0.126
3.688 0.215 0.721 0.421
2.120 0.126 0.421 0.281
study represent a random sample from that population of cracks. The final model for the maximum contrast is
Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
18.079 0.935 3.482 2.034
0.935 0.109 0.199 0.119
3.482 0.199 0.678 0.403
2.034 0.119 0.403 0.273
The model assumptions can be stated succinctly as ε Normal ð0; σ ε Þ and 0 1 00 1 1 β0;b β0 Bβ C BB β C C B 1;b C BB 1 C C B C B C; Σ βb C: B β2;b C Multinormal B @ A @ A 0 @ A β3;b 0
log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ 4. Probability of detection for a population of cracks
β4 log ðvÞlog ðrÞ þ β5 log ðvÞlog ðpÞ þ β6 log ðpÞlog ðrÞ þ β7 log ðvÞ2 þ β8 log ðpÞ2 þβ9 log ðrÞ2 þ ε:
ð2Þ
The random error ε is assumed to have a normal distribution with mean 0 and variance σ ε . For each crack, the coefficients ðβ0;b ; β1;b ; β2;b ; β3;b Þ are assumed to have a multivariate normal distribution independent of ε, with mean μβb ¼ ðβ0 ; β1 ; β2 ; β3 Þ and covariance matrix Σ βb , where μβb is a vector and Σ βb is a 4 4 matrix. These model assumptions can be stated succinctly as ε Normalð0; σ ε Þ and 0 1 00 1 1 β0;b β0 Bβ C BB C C B 1;b C BB β C C B C Multinormal BB 1 C; Σ βb C: B β2;b C BB β 2 C C @ A @@ A A β3;b β3 The parameter estimates, which are calculated by the function lme (in package nlme) in R [17], are given in Table 4 and the matrix elements estimates are given in Table 5.
3.3. Reduced mixed effect model We reduced the full model in (2) based on Akaike's information criterion (AIC) [1] and the t ratios of the regression coefficients. The final model that we used is log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β5 log ðvÞlog ðpÞ þ β9 log ðrÞ2 þ ε:
ð3Þ
We still assumed that ðβ0;b ; β1;b ; β2;b ; β3;b Þ has a multivariate normal distribution independent of ε, but with fixed 0 mean for β2;b and β3;b , i.e., μβb ¼ ðβ0 ; β1 ; 0; 0Þ. The reduced model sets β5, β7 and β8 to be 0. The results for fitting the reduced model are listed in Table 6. The covariance matrix results are given in Table 7.
4.1. Distribution of response variable Under the reduced mixed effect model (3), for given fixed values of v, p and r, the response variable m þ1 has a normal distribution with mean μðv; p; rÞ ¼ β0;b þ β1;b log ðvÞ þ β2;b log ðpÞ þβ3;b log ðrÞ þ β5 log ðvÞlog ðpÞ þ β9 log ðrÞ2
ð4Þ
and variance σ 2 ðv; p; rÞ ¼ X Tl Σ βb X l þ σ 2ϵ :
ð5Þ
Here X l ¼ ð1; log ðvÞ; log ðpÞ; log ðrÞÞT where superscript T indicates vector transpose. Eq. (5) will provide basis for estimating probability of detection (POD) for the population of potential cracks as a function of the experimental factors. 4.2. Calculation of probability of detection Probability of detection (POD) is described in detail in MILHDBK-1823A [14] and Georgiou [8]. In our application, POD is the probability that M from a thermal image sequence exceeds a threshold mT for a crack taken at random from the population of cracks. This probability can be expressed as log ðmT þ1Þ μðv; p; rÞ ð6Þ PODðmT ; v; p; rÞ ¼ PrðM 4mT Þ ¼ 1 Φ σ where ΦðxÞ is the standard normal cumulative distribution function distribution. We can estimate the POD by evaluating Eq. (6) using the estimates of the model parameters in Table 4. The threshold was chosen as mT ¼10. This value was chosen to balance sensitivity with the probability of generating false positives. Fig. 2 displays the output of our mixed-effects model for the POD. The darker the color, the lower the POD value.
232
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Fig. 2. POD plot for pulse length¼ 217. The lighter the color, the higher of the POD.
5. Analysis of the alarm data
Fig. 3. POA plot for pulse length¼ 217. The lighter the color, the larger of the POA.
Some tests ended with an alarm, usually caused by a lock-up of the welder/sample system when the horn can no longer vibrate properly. This is caused by protective feedback loop that terminates power to the welder to prevent damage to the piezoelectric element. Alarm incidents are undesirable because no useful information is obtained from the test. We would like to find test conditions that will also result in a small probability of an alarm. When there was an alarm, the time to alarm was recorded. When there was no alarm, the alarm time was taken to be right censored at the length of the pulse because the unknown alarm time is greater than the pulse length. We fit a log normal regression model to the censored time-to-alarm data to assess the probability of alarm as a function of the test conditions. Let TA denote the time of an alarm. A probability plot of the censored data suggested that TA can be described by a lognormal distribution so that log ðpÞ μa POA ¼ PrðT A rpÞ ¼ Φnor : ð7Þ σa Here μa ¼ βa0 þβa1 v þ βa2 r. Higher order terms were not statistically significant. The restricted maximum likelihood (REML) estimates a a a of the parameters are β^ 0 ¼ 6:027, β^ 1 ¼ 0:0188, β^ 2 ¼ 0:00531, a σ^ ¼ 0:1625. Fig. 3 depicts estimates of POA for pulse length equal to 217 ms.
Fig. 4. Overlay threshold plot for pulse length ¼217. The South-east corner of the plot satisfied both criteria (i.e., POD 4 0:90 and POA o 0:05).
6. Choosing good test conditions
7. Modeling response as a function of effective crack size
A choice of good test conditions needs to fulfill two requirements: a large response for the crack and a small chance of triggering an alarm that would terminate the test and yield no information. Fig. 4 combines the information in the POD plot in Fig. 2 and the POA plot in Fig. 3, showing the test conditions where the probability of detection is expected to be high and probability of alarm is expected to be low. Such overlay threshold plots require specification of acceptable POD and PFA thresholds and we chose those to be POD 4 0:90 and POA o0:05. The south-east region of Fig. 4 meets both conditions, suggesting the test conditions to be used in future tests.
In this section, we will build a statistical model to describe signal amplitude as a function of crack length. This model can then be used to estimate POD as a function of crack length. 7.1. Crack sizes Table 8 lists the crack size information that was obtained by using a method known as acetate replication, a nondestructive technique. Apparent length is the portion of the total wrap-around crack length visible from the same side of the airfoil as the IR camera. Edge length is the portion of the crack length that wraps
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 8 Size information for the 10 cracks.
233
Table 10 Covariance matrix of random effect of full model with size.
Crack
Apparent
Opposite apparent
Edge
Total
Effect
Σ βb
β0
β1
β2
β3
1 2 3 4A 4B 5 6 7A 7B(& C) 7D
0.04 0.1058 0.2903 0.0793 0.2907 0.0935 0 0.2781 0.0962 0.1039
0 0.0304 0.228 0.1875 0 0.1129 0.0718 0.2449 0.1321 0.0503
0.0039 0.021 0.0153 0.0313 0 0.0163 0.0066 0.0156 0.01 0.0126
0.0439 0.1572 0.5336 0.2981 0.2907 0.2227 0.0784 0.5386 0.2383 0.1668
0.0046 0.0226 0.0204 0.0322 0.0054 0.0177 0.0065 0.0205 0.0116 0.0143
β0 β1 β2 β3
21.158 0.929 3.797 2.336
0.929 0.114 0.206 0.128
3.797 0.206 0.705 0.434
2.336 0.128 0.434 0.270
Table 9 Fixed effect estimates of the mixed-effect model with size. βs
Value
Std. error
t-ratio
p-value
β0 β1 β2 β3 β4 β5 β6 β7 β8 β9 β10
3.095 3.517 2.874 1.185 0.662 0.087 0.300 0.213 0.260 0.124 1.122
11.765 1.283 4.713 1.440 0.168 0.087 0.175 0.153 0.476 0.164 0.187
0.263 2.741 0.610 0.823 3.944 0.996 1.716 1.391 0.546 0.755 6.006
0.793 0.007 0.543 0.412 0.000 0.320 0.088 0.165 0.586 0.451 0.000
Table 11 Fixed effect estimates of the reduced mixed-effect model with m as a function of crack size. βs
Value
Std. error
t-ratio
p-value
β0 β1 β5 β6 β9 β10
3.010 1.774 0.637 0.164 0.253 1.127
1.145 0.564 0.108 0.107 0.079 0.186
2.628 3.144 5.874 1.531 3.192 6.048
0.009 0.002 0.000 0.127 0.002 0.000
β3;b (i.e., μβb ¼ ðβ0 ; β1 ; 0; 0Þ). The results for fitting the reduced model are given in Table 11 and the covariance matrix results are given in Table 12.
8. POD as a function of effective crack size 8.1. POD versus crack size
around the edge of the airfoil but is not visible from either side of airfoil. Opposite apparent length is the portion of the crack that would be visible from the opposite side of the airfoil. Especially because there were zeros for some of these reported sizes, we use a convex combination of these recorded sizes to define an effective crack length. We then treat this variable as a fixed effect in the regression model. 7.2. Full mixed effect model with crack size information The full mixed effects model that includes crack length, shown in (8), is similar to model (2). log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β4 log ðvÞlog ðpÞ þ β5 log ðvÞlog ðrÞ þ β6 log ðpÞlog ðrÞ þ β7 log ðvÞ2 þ β8 log ðpÞ2 þβ9 log ðrÞ2 þ β10 log ðLN Þ þ ε:
ð8Þ
Here LN ¼ αLa þ ð1 αÞLe is the effective crack length and α¼0.0185 was chosen to maximize the likelihood. That is, a sequence of models was fit using different values of α between 0 and 1. The value of α that gave the largest value of the likelihood was α¼ 0.0185. The fixed effect parameter estimates of model (8) are given in Table 9. The corresponding estimates of the variances and covariance of the random effect terms are given in Table 10.
It is necessary to estimate POD as a function of crack size for purposes such as planning periodic inspections. In this section, we are going to use model (9) to recalculate the POD function. In contrast to Eq. (6), POD will now also be a function of effective crack size and can be expressed as log ðmT þ 1Þ μðv; p; r; lÞ PODðm; v; p; r; lÞ ¼ Prðm 4 mT Þ ¼ 1 Φ : σ ð10Þ An estimate of POD as a function of effective crack size can be computed by substituting the model estimates from Tables 11 and 12 into model (9). Fig. 5 shows the POD estimate as a function of effective crack length for different test conditions (Vibration Amplitude, Pulse Length and Trigger Force, respectively). These particular conditions were chosen to provide examples of different POD values. The inside tic marks at the bottom of the plot show the measured effective sizes of the cracks in the data set. The POD estimates to the right of effective size 0.03 involve extrapolation. Both of the full model estimates and reduced model estimates are shown in the plot indicating little difference between the two models when estimating POD. In the following sections wee use the reduced model for its simplicity. 8.2. Bootstrap confidence interval of POD estimation
7.3. Reduced mixed effect model with crack size information Similar to what we did in Section 3.3, the full model can be reduced to a smaller model based on AIC. The model suggested by this criterion is log ðM þ 1Þ ¼ β0;b þ β1;b log ðvÞ þβ2;b log ðpÞ þ β3;b log ðrÞ þ β5 log ðvÞlog ðrÞ þ β6 log ðpÞlog ðrÞ þβ9 log ðrÞ2 þβ10 log ðLN Þ þ ε:
ð9Þ
We again assume that ðβ0;b ; β1;b ; β2;b ; β3;b Þ has a multivariate normal distribution independent of ε, but with fixed 0 mean for β2;b and
The complicated model fitting procedure and calculation of POD make the direct computation of a confidence interval of POD difficult. Bootstrap methods provide a convenient alternative for such complicated procedures as described in Efron and Tibshirani [5]. The bootstrap procedures that we used are fully parametric. That is, we simulate bootstrap samples from our assumed parametric model and the REML estimates [3]. Fig. 6 shows the result for both the percentile and the BCA (bias-corrected accelerated) pointwise nonparametric bootstrap confidence intervals for POD
234
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
Table 12 Estimated covariance matrix of random effect for the reduced model with m as a function of size. Σ βb
β0
β1
β2
β3
β0 β1 β2 β3
20.302 0.884 3.633 2.257
0.884 0.109 0.197 0.123
3.633 0.197 0.674 0.419
2.257 0.123 0.419 0.263
1.00
POD vs. Effective Crack Length Contrast Threshold = 10
Coverage Probability
0.95
● ●
● ●
0.90
●
●
●
●
●
●
●
●
0.85
● ● ●
0.80
●
●
●
●
BCA 95% Confidence Interval, 10 Samples Percentile 95% Confidence Interval, 10 Samples BCA 95% Confidence Interval, 20 Samples Percentile 95% Confidence Interval, 20 Samples BCA 95% Confidence Interval, 40 Samples Percentile 95% Confidence Interval, 40 Samples
0.75
0.8
0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040
Partial Model Full Model 0.6
● ●
1.0 (60,217,60)
●
Crack Length Fig. 7. Coverage probability of POD Confidence Intervals calculated by the BCA and the percentile bootstrap method for 10, 20 and 40 crack samples. The solid line is by the BCA method and the dotted line is by the percentile method. For the range of the crack size of our sample, the BCA method has better coverage than the percentile method.
POD
(35,150,35)
0.4
versus effective crack length in test condition (35, 150, 35). 5000 bootstrap samples were used in the calculation.
0.2
8.3. Coverage Probabilities of the Percentile and the BCA Confidence Interval
(10,83,10) 0.0 0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
Effective Crack Length
Fig. 5. POD as a function of effective crack length. The solid line is for the reduced model and the dotted line is for the full model. The three test condition parameters (Vibration Amplitude, Pulse Length and Trigger Force) are indicated in the parentheses near the lines.
1.0
Although bootstrap procedures are expected to provide actual coverage probabilities that are closer to nominal than simpler methods such as the Wald's approximate method, the actual coverage probability, in general, will only be approximate. We used simulations to compare the coverage probabilities of the two bootstrap confidence interval calculation methods: the sample percentile method and the more refined the BCA method using 2000 new data sets were simulated from the fitting result of the reduced model and the percentile and the BCA bootstrap methods were used to calculate confidence intervals. The observed coverage probability is defined as Coverage ¼
0.8
0.6 POD
BCA Confidence Interval Percentile Confidence Interval 0.4
0.2
0.0 0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
Effective Crack Length
Fig. 6. 95% confidence intervals based on the percentile and the BCA bootstrap methods for test condition (Vibration Amplitude ¼35, Pulse Length¼ 150 and Trigger Force¼35). The solid line in the middle is the POD estimate from the mixed effect model. The dotted lines are the pointwise confidence intervals, where the large dotted line is for the BCA method and the small dotted line is for the percentile method.
NI N
ð11Þ
where NI denotes the number of confidence intervals that include the true value of POD (which is calculated in the model from which the data sets are simulated) and N is the total number simulated data sets. The simulation results of the coverage probability are shown in Fig. 7, for data sets with 10, 20 and 40 experimental units with three different symbols. The solid lines are for the BCA method and the dotted lines are for the percentile method. The BCA shows closer to nominal converge probabilities than the percentile method, as expected. The coverage probability of both the BCA method and the percentile method become closer to nominal when the sample size increases. The BCA method with 40 experimental units is the closest to nominal among all six evaluations. When sample size is small, (i.e., 10), the coverage probabilities are far away from nominal for larger cracks. As the sample size increases, the difference in coverage between small cracks and large cracks becomes smaller.
9. Conclusion and future work This paper has shown how to find good vibrothermography inspection test conditions for detecting cracks in fan blades and
C. Gao et al. / Reliability Engineering and System Safety 131 (2014) 229–235
how to estimate POD as a function of the effective size of a crack. The same general approach could be use for other kinds of inspections where tuning parameters need to be set. There are a number of areas for possible future research in this area.
It would be useful to develop methods for quantifying and displaying uncertainty in the overlay threshold plots like Fig. 2.
In our work, the standard maximum contrast response was
used. It would be useful to explore the use of other response variables in an attempt to use more completely the information in the sequence of image produced by a vibrothermography test system. For example, Li et al. [11] used a signal-to-noise ratio criterion after a matched filter as a response from a different vibrothermography inspection experiment. Gao and Meeker [7] use a method based on principal components analysis of a set of coefficients from a series of fitted regression models to extract a scalar measure of increase in temperature from the sequence-or-image data. With the use of diffuse prior information, Bayesian methods will provide results that are similar to the REML method used in our study. The Bayesian method provides an attractive method to construct confidence intervals for complicated functions of the parameters (such as POD) and results in other areas of application have shown that intervals computed in this manner have good frequentist coverage properties. Li et al. [12] use Bayesian methods with a mixed effect model to estimate POD based on ultrasonic inspection of synthetic hard alpha inclusions in a disk forging.
Acknowledgments This material used in the paper is based upon work supported by the Air Force Research Laboratory under Contract # FA8650-04C-5228 at Iowa State University's Center for NDE.
235
References [1] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control 1974;19(6):716–23. [2] Badia FG, Berrade MD, Campos CA. Optimal inspection and preventive maintenance of units with revealed and unrevealed failures. Reliab Eng Syst Saf 2002;78:157–63. [3] Bartlett MS. Properties of sufficiency and statistical tests. Proc R Soc Lond Ser A Math Phys Sci 1937;160(901):268–82. [4] Cronvall O, Simola K, Männistö I, Gunnars J, Alverlind L, Dillström P, et al. A study on the effect of flaw detection probability assumptions on risk reduction achieved by non-destructive inspection. Reliab Eng Syst Saf 2012;105:90–6. [5] Efron B, Tibshirani RJ. An introduction to the bootstrap. Chapman & Hall/CRC; 1994. [6] Favro LD, Han X, Ouyang Z, Sun G, Sui H, Thomas RL. Infrared imaging of defects heated by a sonic pulse. Rev Sci Instr 2000;71(6):2418–21. [7] Gao C, Meeker WQ. A statistical method for crack detection using vibrothermography inspection data. Qual Technol Quant Manag 2012;9:59–77. [8] Georgiou GA. Probability of detection (pod) curves derivation, applications and limitations. Technical report, Jacobi Consulting Ltd, 2006. Downloaded from 〈http://www.hse.gov.uk/research/rrpdf/rr454.pdf〉 on May 6; 2011. [9] Hellier C. Handbook of nondestructive evaluation. New York: McGraw-Hill Professional; 2001. [10] Henneke EG, Jones TS. Detection of damage in composite materials by vibrothermography. ASTM Spec Tech Publ 1979;1(696):83–95. [11] Li M, Holland SD, Meeker WQ. Statistical methods for automatic crack detection based on vibrothermography sequence-of-images data (with discussion). Appl Stoch Models Bus Ind 2010;26:481–95. [12] Li M, Meeker WQ, Thompson RB. Physical model-assisted probability of detection of flaws in titanium forgings using ultrasonic nondestructive evaluation. Technometrics 2014;56:78–91. [13] Mayton D, Spencer F, Alvarez C. Characterizing the effects of sonic IR variables on turbine disk inspection using a design of experiments approach. Rev Progr Quant Nondestruct Eval 2005;24:642–9. [14] MIL-HDBK-1823A. Nondestructive evaluation system reliability assessment, Standardization Order Desk, Building 4D, 700 Roberts Avenue, Philadelphia, PA 19111-5094, 2009. Downloaded from 〈http://www.statisticalengineering. com/mh1823/MIL-HDBK-1823A%282009%29.pdf〉 on May 6; 2011. [15] Mix PE. Introduction to nondestructive testing: a training guide. New York: Wiley-Interscience; 2005. [16] Pinheiro JC, Bates DM. Mixed-effects models in S and S-Plus. New York: Springer; 2000. [17] R Development Core Team, R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2008. [18] Shull PJ. Nondestructive evaluation: theory, techniques, and applications. New York: Marcel Dekker, Inc.; 2002. [19] Simola K, Pulkkinen U. Models for non-destructive inspection data. Reliab Eng Syst Saf 1998;60:1–12. [20] Venables WN, Ripley BD. Modern applied statistics with S. Fourth edition. New York: Springer; 2002. [21] Xiong JJ, Shenoi RA. A durability model incorporating safe life methodology and damage tolerance approach to assess first inspection and maintenance period for structures. Reliab Eng Syst Saf 2009;94:1251–8.