- Email: [email protected]
A competing risk model for the reliability of cylinder liners in marine Diesel engines
ARTICLE IN PRESS Reliability Engineering and System Safety 94 (2009) 1299–1307
Contents lists available at ScienceDirect
Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress
A competing risk model for the reliability of cylinder liners in marine Diesel engines D. Bocchetti a, M. Giorgio b, M. Guida c, G. Pulcini d, a
Grimaldi Group, Naples, Italy Department of Aerospace and Mechanical Engineering, Second University of Naples, Aversa (NA), Italy Department of Information Engineering and Electrical Engineering, University of Salerno, Fisciano (SA), Italy d Istituto Motori, National Research Council—CNR, Naples, Italy b c
a r t i c l e in fo
abstract
Article history: Received 13 November 2007 Received in revised form 29 January 2009 Accepted 29 January 2009 Available online 10 February 2009
In this paper, a competing risk model is proposed to describe the reliability of the cylinder liners of a marine Diesel engine. Cylinder liners presents two dominant failure modes: wear degradation and thermal cracking. The wear process is described through a stochastic process, whereas the failure time due to the thermal cracking is described by the Weibull distribution. The use of the proposed model allows performing goodness-of-fit test and parameters estimation on the basis of both wear and failure data. Moreover, it enables reliability estimates of the state of the liners to be obtained and the hierarchy of the failure mechanisms to be determined for any given age and wear level of the liner. The model has been applied to a real data set: 33 cylinder liners of Sulzer RTA 58 engines, which equip twin ships of the Grimaldi Group. Estimates of the liner reliability and of other quantities of interest under the competing risk model are obtained, as well as the conditional failure probability and mean residual lifetime, given the survival age and the accumulated wear. Furthermore, the model has been used to estimate the probability that a liner fails due to one of the failure modes when both of these modes act. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Competing risk model Reliability Diesel engine cylinder liners Wear process Thermal cracking
1. Introduction Modern marine Diesel engines must satisfy stringent reliability and availability requirements. Experience shows that the more critical components of heavy-duty Diesel engines are the cylinder liners. These components are very costly and their in-service failures can lead to heavy economic losses. As such, estimating their reliability is crucial when performing long-term economic evaluations and to plan maintenance activities. Unfortunately, reliability estimation of cylinder liners can be difficult to perform due to the scarcity of failure data, and large experimental costs prevent the running of ad hoc reliability tests. As a consequence and for economic reasons, the liners undergo extensive maintenance, aimed at minimizing the number of inservice failures. As such, it is important to make use of models and methodologies that can allow good reliability estimates to be performed based on the available data and all the other technological and experimental information usually possessed by the analyst. Cylinder liners present two dominant failure modes: wear degradation and thermal cracking. In high power marine Diesel
Corresponding author. Tel.: +39 0817177113; fax: +39 081 2396097.
E-mail address: [email protected] (G. Pulcini). 0951-8320/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2009.01.010
propulsion engines, wear usually occurs in the top region of the inner surface of the cylinder liner, where the maximum mechanical and thermal load occurs. Many studies of Diesel engine wear agree that the dominant wear mechanism in this region is due to the high quantity of abrasive particles on the piston surface, generated by the combustion of heavy fuels and oil degradation (soot) (see Refs. [1–3]). The abrasive wear mechanism can be considered as a three-body contact mode. If the lubricant film thickness is less than the soot particle size, then the soot can cause an abrasive action with the liner metal surfaces. Indeed, soot particles are harder than the corresponding engine parts, which further increase the process of abrasion. In addition to abrasive wear, corrosive wear (due to sulphuric acid, nitrous/nitric acids, and water) can also be observed. The thermal cracking of a liner is due to a fatigue mechanism originated by repeated thermal shocks. A thermal shock is caused by a rapid change in the temperature of the cooling fluid that laps the external surface of the liner, and is often aggravated by scale and corrosion in the cooling water spaces as a result of inadequate chemical treatment of the cooling water. These temperature changes occur mainly during the vessel manoeuvre when it is difficult to maintain constant the cooling fluid temperature. As a consequence of the repeated thermal shocks, a small crack can arise in the external surface of the liner and, with time, it propagates towards the inside.
ARTICLE IN PRESS 1300
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
In this paper, a competing risk model is proposed to describe the failure process of a cylinder liner under the hypotheses that (a) the liner fails when the first of the competing failure mechanisms reaches a failure state, (b) the failure mechanisms proceed independently of each other, and (c) each failure mode has a known life distribution. Based on wear measurements taken from a large set of liners at different ages, the wear failure mode is treated as a degradation process for which the failure occurs when a fixed warranty limit is crossed. On the contrary, in the absence of any information on the underlying process of thermal shocks leading to fatigue cracking, the thermal crack mode is treated as a catastrophic failure mode, and the time to failure due to the thermal shocks are modeled via a Weibull distribution. The proposed model has been applied to a real data set referring to the cylinder liners of Sulzer RTA 58 engines, which equip twin ships of the Grimaldi Group. Estimates of the liner reliability and failure rate are obtained based on the competing risk model approach. The probability of failure of the liner is caused by one of the failure modes, when both of these modes act. By its very nature, using the proposed model allows the analyst to: (a) justify the choice of the model on the basis of strong experimental and technological arguments; (b) estimate the reliability of liners using both failure and wear data; (c) use information collected during recurrent maintenance inspections to improve the accuracy of reliability estimates and to perform estimates conditioned to the state of the liners; and (d) recognize the failure modes hierarchy for any given age and wear level of the liner, thus providing the analyst with information about the criticality of the failure mechanisms.
Fig. 1. The cross-section of the cylinder of SULZER RTA 58 engine.
tions are survival times. In three cases (i ¼ 12, 25, and 28), the liner was replaced due to a large accumulated wear, and hence the survival time ti coincides with the time ti;ni of the last inspection. The total cumulated operating time of all the 33 liners was of 913,894 h.
3. Model assumptions 2. The case study
3.1. Independence of failure modes
The failure data set consists of m ¼ 33 cylinder liners of 8cylinder SULZER RTA 58 engines [4], equipping twin ships of the Grimaldi Group. These data were collected from the stage each liner was initially installed, between the dates of January 1, 1994 and February 28, 2007. The liners were recurrently inspected and the accumulated wear was measured. Measurements were taken by positioning a caliper inside the cylinder bore at the position of the Top Dead Center (TDC) (Fig. 1): the point of maximum wear [2]. The sensitivity of the caliper was 0.05 mm, and as such each measure was rounded up to the nearest 0.05. Furthermore, when a liner failed due to thermal cracking, it was replaced with a new one. Thus, the wear data set consisted of the number ni (i ¼ 1, y, 33) of inspections performed on each liner, the age ti,j (j ¼ 1, y, ni) (in operating hours) of the liner at the time of inspection j, and the wear Wi,jW(ti,j) (in mm) measured during P each inspection (see Table 1). A total of N ¼ m i¼1 ni ¼ 64 measures of wear were recorded. The liner i ¼ 32 failed due to thermal cracking before the first inspection was performed, and thus no wear data was subsequently recorded from it. Thermal shocks were not recorded, even if the cooling fluid temperature was monitored, such that the data concerning the thermal crack (given in Table 1) consist only of the age ti (i ¼ 1, y, 33) at which either the failure occurred (di ¼ 1) or did not occur (di ¼ 0), due to thermal shock. A total of 22 failures by thermal cracking were observed, and the remaining 11 observa-
The two failures mechanisms, namely the wear and the thermal crack, proceed independently of each other (at least until a failure occurs). This assumption is justified by the fact that the physical mechanisms that cause or influence the evolution of each failure mode are quite different and unrelated. In fact, wear that occurs on the inner surface of the liner cannot affect the mechanical strength of the liner. Indeed, the maximum admissible wear (4 mm) is negligible with respect to the initial liner thickness (100 mm). Furthermore, the point at which maximum wear occurs is far enough from the point at which thermal crack initiation occurs and subsequently propagates towards the inside. Also, the changes in the temperature of the cooling fluid, which are the driving mechanism for thermal cracking, do not influence the wear formation process, the size of abrasive particles, or lubrication oil degradation. The independence assumption implies that the probability of a liner failure in (t,t+dt) from one of the two failure mechanisms, given survival to time t, is the same as to whether that mechanism were the only failure mechanism acting on the liner. 3.2. The wear model Using the wear measurements taken on the liners at different ages, the wear failure mode is treated as a degradation process for
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
1301
Table 1 Wear and thermal crack data of cylinder liners. i
ni
ti,1 (h)
Wi,1 (mm)
ti,2 (h)
Wi,2 (mm)
ti,3 (h)
Wi,3 (mm)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
3 2 2 2 3 3 3 3 3 4 3 4 3 2 1 1 1 1 1 1 1 1 1 1 4 1 1 3 2 1 1 0 1
11,300 11,360 11,300 12,300 14,810 9700 10,000 6860 2040 7540 8510 18,320 10,000 9350 13,200 7700 7700 8250 3900 7700 9096 19,800 10,450 12,100 12,000 8800 2200 33,000 8800 8250 18,755
0.90 0.80 1.50 1.00 1.90 1.10 1.20 0.50 0.40 0.50 0.80 2.20 2.10 0.85 2.00 1.05 1.60 0.90 1.15 1.20 0.50 1.60 0.40 1.00 1.95 1.40 0.40 2.90 0.50 0.70 1.15
14,680 17,200 21,970 16,300 18,700 19,710 30,450 17,200 12,580 8840 14,930 25,310 16,620 15,970
1.30 1.35 2.00 1.35 2.25 2.60 2.75 1.45 2.00 1.10 1.45 3.00 2.75 1.20
31,270
2.85
28,000 30,450 37,310 24,710 16,620 9770 21,560 37,310 30,000
2.75 3.00 3.05 2.15 2.35 1.15 1.90 3.70 3.60
27,300
2.70
49,500
3.15
38,500 27,500
3.25 2.15
55,460
4.10
8490
0.95
which the degradation warranty limit Wmax ¼ 4 mm is considered [5]. The wear process is then modeled through a point stochastic process, by assuming (as suggested in Ref. [6]) that: (a) the wear of a cylinder liner at a given age t is the accumulation of isolated damage (elementary wear) of equal magnitude c that has occurred by age t; (b) the events that cause an elementary wear occur randomly in time; and (c) the probability that damage occurs in a given time interval is independent of the previously accumulated damage.
0
0
a0 40; 0ob o1.
However, since the wear process cannot grow unlimitedly, a stochastic process that bounds the expected number of wearing events offers a more realistic modeling approach. Considering this, the wearing process can now be described through a nonhomogeneous Poisson process (NHPP) with intensity function:
lðtÞ ¼ a expðbtÞ;
a40; 1obo1
(1)
often referred to as log-linear process (LLP) [8]. Based on the LLP assumption, the expected number of wear events up to the operating time t, say MðtÞ E½NðtÞ, is given by MðtÞ E½NðtÞ ¼ ða=bÞ½expðbtÞ 1
(2)
Note that, when bo0, M(t) tends to a/b as the age t of the liner increases. Assumptions are made in terms of an independent increments process and of constant elementary damages in order to avoid
Wi,4 (mm)
16,300
2.10
45,000
3.95
56,120
4.05
ti (h)
di
36,370 28,930 27,970 21,830 39,450 39,500 42,320 25,200 27,750 25,680 29,900 45,000 37,860 25,100 18,270 18,650 26,770 16,870 11,600 14,300 14,596 31,900 25,300 23,720 56,120 20,950 16,738 55,460 28,100 8250 31,330 5430 16,790
1 1 1 1 0 1 0 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 0 0 0 1 0 1 0 1 1 1
that some estimation issues, such as those discussed later in Section 4.1, become unsolvable. However, despite these simplifications, the proposed model proved to be suitable for this case study. On the basis of the above assumptions, the wear accumulated at the operating time t is modeled by the stochastic process: WðtÞ ¼ cNðtÞ
(3)
so that the probability distribution of the wear w(t1, t2) accumulated in the time interval (t1, t2) is P½wðt1 ; t 2 Þ ¼ w ¼ P½Nðt 1 ; t 2 Þ ¼ w=c
Ref. [6] describes the wearing events process N(t) through a power law process (PLP) [7] with decreasing intensity function:
l0 ðtÞ ¼ ðb0 =a0 Þðt=a0 Þb 1 ;
ti,4 (h)
¼
½Mðt 2 Þ Mðt 1 Þw=c exp½Mðt 2 Þ þ Mðt 1 Þ ðw=cÞ!
(4)
where N(t1, t2) denotes the number of wear events in the interval (t1, t2) and w ¼ 0, c, 2c,y . For each NHPP, Var[N(t)] ¼ M(t). Therefore, the mean and variance of the wear accumulated until t are, respectively E½WðtÞ ¼ cMðtÞ and Var½WðtÞ ¼ c2 MðtÞ
(5)
Since the time at which the wear level crosses the warranty limit Wmax ¼ 4 mm is considered as a failure time, the wear mode reliability of a liner, defined as the hypothetical probability that, in absence of the other failure modes, the wear accumulated until t does not exceed the 4 mm limit [9], is given by RW ðtÞ ¼ P½WðtÞpW max ¼ P½NðtÞpK max ¼
K max X
½MðtÞk exp½MðtÞ ¼ 1 IG½MðtÞ; K max þ 1 k! k¼0
(6)
where K max ¼ int½W max =c is the greatest integer equal to or less Ry than Wmax/c, and IGðy; kÞ ¼ 0 xk1 expðxÞdx=GðkÞ is the incom-
ARTICLE IN PRESS 1302
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
plete Gamma function. The probability density function (pdf) is f W ðtÞ ¼
a½MðtÞ
GðK max þ 1Þ
exp½bt MðtÞ;
t40
(7)
RW ðtjWðt 1 ÞÞ ¼ P½WðtÞpW max jWðt 1 Þ ¼ P½Nðt 1 ; tÞpDK max DX K max k¼0
½MðtÞ Mðt 1 Þk exp½MðtÞ þ Mðt 1 Þ k!
¼ 1 IG½MðtÞ Mðt 1 Þ; DK max þ 1
(8)
where DK max ¼ int½fW max Wðt 1 Þg=c is the greatest integer equal to or less than {WmaxW(t1)}/c. Fig. 2 shows the measured wear accumulated by the cylinder liners, where the measured points for each liner have been linearly interpolated for graphical representation. The wear behavior is in keeping with the literature results (see, for example, Ref. [10]) and the technological characteristics: the early phase is a wear accommodation period in which the rate of wear decreases quickly, followed by a phase in which the rate decreases more slowly. In the present application, the accommodation phase occupies the large part of the entire observation time, since the warranty clauses impose that a cylinder liner is replaced before it accumulates a relatively ‘‘low’’ wear level (4 mm), in order to avoid potentially catastrophic and very expensive failures. Moreover, the observed (interpolated) processes are closely interwoven, thus indicating that the rate of wear of an individual liner changes continually, and that these changes are of a random nature. 3.3. The thermal crack model In absence of information regarding the underlying process of thermal shocks, leading to fatigue cracking, the thermal crack mode is treated as a hard (or catastrophic) failure mode, and its reliability at the age t (i.e. the hypothetical probability that a thermal crack failure occurs after t in absence of the wear mode) is assumed to be Weibull RTC ðtÞ ¼ exp½ðt=aÞb ;
a; b40
(9)
This assumption is justified on the basis of theoretical considerations since fatigue life data is often shown to be adequately analyzed using the Weibull distribution (see Refs. [11–13]), and is supported by graphical analysis. In particular, the graphical analysis was performed by plotting, on Weibull paper, the Kaplan–Meier estimate [14] of the reliability
Measured wear [mm]
4.0 3.5 3.0 2.5 1 4 7 10 13 16 19 22 25 28 31
2.0 1.5 1.0 0.5
2 5 8 11 14 17 20 23 26 29
3 6 9 12 15 18 21 24 27 30 33
0.0 0
R¯ i R¯ TC ðti Þ ¼
Y rl 1
10,000
20,000
30,000
40,000
Operating time [h] Fig. 2. The observed wear processes.
50,000
60,000
(10)
rl
l2Si
Given the wear W(t1)oWmax accumulated by time t1, the conditional probability that the warranty limit is not exceeded before t4t1 is given by
¼
at any catastrophic failure time ti (i:di ¼ 1) versus ti:
K max
where Si is the set of liners whose catastrophic failure occurred at tlpti, and rl is the reverse rank, i.e., the number of failure or censoring times no less than rl (see Table 2, where the catastrophic failure and censoring times were ordered for convenience reason). For example, S4 ¼ (4, 15, 16, 18, 19, 20, 21, 27, 32, 33) and r4 ¼ 23. As depicted in the Weibull plot of Fig. 3, the points roughly follow a straight line and give no obvious evidence that the thermal crack data do not fit a Weibull distribution. Furthermore, the high determination coefficient, say R2 ¼ 0.95, shows a high goodness of fit. Table 2 also gives both the estimated standard deviation of the estimated reliability [15], based on the following function:
sfR¯ i g ¼ R¯ i
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ½r l ðr l 1Þ1
(11)
l2Si
and the approximate g ¼ 0.80 confidence interval on RTC(ti), which is based on logit transformation [16] and is given by
R¯ i R¯ i ; R¯ i þ ð1 R¯ i Þs R¯ i þ ð1 R¯ i Þ=s
(12)
where s ¼ expfz1g=2 sfR¯ i g=½R¯ i ð1 R¯ i Þg and z1g/2 is the (1g/2) quantile of the Normal standard distribution. As shown in Fig. 3, the fitted line falls within the 80% confidence intervals. Note that the endpoints of Eq. (12) always lie between 0 and 1. Table 2 Kaplan–Meier estimate and 80% confidence interval of reliability for thermal crack mode. Liner i
ti
di
ri
R¯ i
sfR¯ i g
80% confidence interval
32 19 20 21 27 33 18 15 16 26 4 24 14 8 23 10 17 9 3 30 29 2 11 31 22 1 13 5 6 7 12 28 25
5430 11,600 14,300 14,596 16,738 16,790 16,870 18,270 18,650 20,950 21,830 23,720 25,100 25,200 25,300 25,680 26,770 27,750 27,970 28,050 28,100 28,930 29,990 31,330 31,900 36,370 37,860 39,450 39,500 42,320 45,000 55,460 56,120
1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 1 0 1 1 0 1 1 1 1 1 1 0 0 1 0 0 0 0
33 32 31 30 29 28 27 26 25
0.9697 0.9394 0.9091 0.8788 0.8485 0.8182 0.7879 0.7576 0.7273
0.0298 0.0298 0.0415 0.0500 0.0567 0.0622 0.0669 0.0709 0.0742
(0.8970, 0.9916) (0.8879, 0.9681) (0.8401, 0.9501) (0.7990, 0.9297) (0.7609, 0.9079) (0.7247, 0.8850) (0.6898, 0.8612) (0.6559, 0.8367) (0.6228, 0.8116)
23
0.6957
0.0774
(0.5885, 0.7851)
20 19 18
0.6609 0.6261 0.5913
0.0810 0.0839 0.0861
(0.5508, 0.7559) (0.5141, 0.7260) (0.4782, 0.6955)
16 15
0.5543 0.5174
0.0883 0.0898
(0.4403, 0.6629) (0.4034, 0.6296)
13 12 11 10 9 8
0.4776 0.4378 0.3980 0.3582 0.3184 0.2786
0.0913 0.0920 0.0918 0.0909 0.0890 0.0864
(0.3638, 0.5937) (0.3254, 0.5569) (0.2880, 0.5193) (0.2517, 0.4808) (0.2164, 0.4414) (0.1821, 0.4011)
5
0.2229
0.0852
(0.1325, 0.3501)
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
under the competing risk model, is
0.80 0.70 0.60 0.50 0.40
hC ðtÞ ¼ hW ðtÞ þ hTC ðtÞ ¼ b1
þ
0.30
a½MðtÞkmax exp½bt MðtÞ Gðkmax þ 1Þf1 IG½MðtÞ; kmax þ 1g
b t
(17)
a a
The mean lifetime can be obtained by numerical integration of the liner reliability (13): Z 1 RC ðtÞdt (18) EC ½T ¼
0.20 1 - RTC (t)
1303
0
0.10
Given that the thermal crack has not occurred until t1 and the accumulated wear at t1 is Wðt 1 ÞoW max , the conditional probability that the liner does not fail before t4t1 is given by " b # t1 b t RC ðtjt 1 ; Wðt 1 ÞÞ ¼ f1 IG½MðtÞ; Dkmax þ 1gexp
0.05
a
0.02
00 ,0
00 ,0
The mean residual life, given the survival age t1 and the accumulated wear W(t1), is given by the integration of (19) over the interval (t1,N). Furthermore, the conditional pdfs at tXt1 for the wear and thermal crack mode are, respectively
50
Operating time (in log scale)
f C;W ðtjt 1 ; Wðt 1 ÞÞ ¼
Fig. 3. Graphical analysis of thermal crack data on the Weibull plot.
a½MðtÞ Mðt 1 ÞDkmax GðDkmax þ 1Þ " exp bt MðtÞ þ Mðt 1 Þ
3.4. The competing risk model
f C;TC ðtjt 1 ; Wðt 1 ÞÞ ¼
Under the assumption of independence, the reliability, under the competing risk model, is the product of the failure mode reliability functions (6) and (9) [9,17,18]: RC ðtÞ ¼ RW ðtÞRTC ðtÞ
(13)
The probability density functions of each failure mode in the presence of both the modes are given by f C;W ðtÞ ¼ f W ðtÞRTC ðtÞ ¼
a½MðtÞkmax
Gðkmax þ 1Þ b t b1
exp½bt MðtÞ ðt=aÞb
f C;TC ðtÞ ¼ f TC ðtÞRW ðtÞ ¼ a a " # b t f1 IG½MðtÞ; kmax þ 1g exp
a
f C ðtÞ ¼ f C;W ðtÞ þ f C;TC ðtÞ
(14)
(15)
Thus, from a computational viewpoint, it is more convenient to obtain first the failure probability FC,W(t) by integrating numerically its pdf fC,W(t), and then derive FC,TC(t) from F C;TC ðtÞ ¼ 1 RW ðtÞRTC ðtÞ F C;W ðtÞ
and b t b1
a a
b t
a
þ
b # t1
a
" b # t1 b t exp
a
a
f1 IG½MðtÞ Mðt 1 Þ; Dkmax þ 1g
(20)
The conditional failure probabilities F C;W ðtjt 1 ; Wðt 1 ÞÞ and F C;TC ðtjt 1 ; Wðt 1 ÞÞ are then given by the (numerical) integration of (20). Their asymptotic value, for t-N, provides the probability that wear degradation and thermal cracking failure occur, respectively, given the state of the liner at time t1.
4. Estimation procedure
The associated failure probabilities, say FC,W(t) and FC,TC(t), are given by the (numerical) integration of (14), and their asymptotic value, for t-N, is less than 1 and provides the proportion of failures caused by each failure mode during the whole life of a liner. The liner pdf and cumulative distribution function (Cdf), under the competing risk model, are thus given by the following, respectively:
F C ðtÞ ¼ F C;W ðtÞ þ F C;TC ðtÞ
a
(19)
40
30
,0
00 ,0 20
10
5,
,0
00
00
0
0.01
00
Approx. 80% Confidence Interval Non Parametric Estimate Weibull Fit
(16)
Under the independence assumption, the failure mode hazard functions hC,W(t) and hC,TC(t) equal the marginal hazard functions hW(t) and hTC(t), respectively, so that the liner hazard function,
Under the assumption of Section 3, the likelihood function relative to the observed data is 9 8 ni m Y
i¼1 m b X ti i¼1
#
a
i¼1
a a
(21)
where wi;j ¼ Wðt i;j Þ Wðt i;j1 Þ is the wear accumulated by the liner i during the inspection interval (ti,j1, ti,j) (when W(ti,0)0 and ti,00), and T i t i;ni is the age of the liner i at the time of the last inspection. Since the stochastic models for the two failure modes are indexed by uncommon parameters, then the parameter estimation can be performed separately for each mode. 4.1. Estimate of wear model parameters For the wear model (3), both the parameters a and c of the intensity function (1) act as multiplicative factors for the accumulated wear. Hence, the maximum likelihood estimation (MLE) procedure cannot be used to estimate a and c individually.
ARTICLE IN PRESS 1304
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
Thus, to overcome this estimation problem, the following fivestep estimation procedure, that combines MLE and first moments estimation methods, is used (see Appendix A for the mathematical proof): 1. Obtain the MLE b^ by solving numerically the equation
(22)
3. Estimate the expected wear accumulated by each liner i during each interval (ti,j1, ti,j) i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni (24) and compute the N ¼ y^ k ¼
^ i;j wi;j E½w qffiffiffiffiffiffiffiffiffiffiffiffiffi ; ^ i;j E½w
Pm
i¼1 ni
values of
i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni ; k ¼ j þ
i1 X
nl
l¼0
(25) 4. Estimate the elementary wear c as PN ðy^ Ave½y^ k Þ2 c^ ¼ Var½y^ k ¼ k¼1 k N1 P ^ ^ where Ave½y^ k ¼ N k¼1 yk =N is the sample mean of yk . 5. Estimate the LLP parameter a as ^ c^ a^ ¼ y^ b=
(26)
(27)
Note that, since E½WðtÞ ¼ y½expðbtÞ 1, its estimate is a MLE. 4.2. Estimate of thermal crack model parameters For the thermal crack model (9), the MLE procedure for Weibull randomly censored sample [15] is used. In particular, the MLE of the parameter b is the numerical solution of Pm b Pm ti lnðti Þ 1 d lnðti Þ þ i¼1 i¼1 ¼0 (28) Pmi Pm b i¼1 di
t
i¼1 i
and the MLE of a is 2P
a^ ¼
31=b^ b^ m i 5 4 Pi¼1 m i¼1 di
t
c^ ¼ 0:1810 mm
The MLE of b is negative, thus indicating that the rate of wear decreases with the liner age. The sample mean of y^ k is equal to 0.0351. This value is much less than the sample standard deviation sðy^ k Þ ¼ 0:4254, thus giving confidence that the estimate of c is sufficiently accurate. ^ and c^ have been obtained without The estimates a^ , b, constraining the ratios wi,j/c, which give the number of wear events experienced by the liner i during the inspection interval (ti,j1, ti,j), to be integer values. Such a lack of restriction does not imply any mathematical issues, even though Eq. (A.1) no longer represents a probability statement. However, because the isolated damages are not of equal magnitude and the data set consists of rounded measures of wear, the lack of restriction on the ratios wi,j/ c does not produce misleading estimates. For a comparison, the wear data are analyzed also within the 0 0 0 PLP model with intensity function l ðtÞ ¼ ðb =a0 Þðt=a0 Þb 1 proposed in Ref. [6]. The resulting MLE of the PLP parameters are 0 0 a^ ¼ 776:5 h and b^ ¼ 0:7458, and the estimate of the elementary 0 wear is c^ ¼ 0:1735 mm. The goodness of fit of the wear models is compared through: (a) the estimated log-likelihoods relative to ^ c^ Þ and lða^ 0 ; b^ 0 ; c^ 0 Þ, and (b) the mean squared the wear data, say lða^ ; b; error (MSE), defined as Pm Pni 2 ^ i¼1 j¼1 ðwi;j E½wi;j Þ MSE ¼ m where
The accuracy of the estimate of c, and hence of a, depends ^ i;j and on the sample size N. A larger N on the accuracy of E½w provides a more accurate estimate of c, as given by Eq. (26). If the sample mean of y^ k is close to 0, then we can be confident that the estimate of c is accurate.
b
b^ ¼ 2:423 105 =h;
(30)
where W(Ti) is the accumulated wear of the liner i by the last inspection. 2. Obtain the MLE of y ¼ ca/b from Pm WðT i Þ (23) y^ ¼ Pm i¼1 ^ i¼1 ½expðbT i Þ 1
^ Þ expðbt ^ ^ i;j ¼ y^ ½expðbt E½w i;j i;j1 Þ;
On the basis of the estimation procedure of Section 4.1, the estimate of the wear process parameters relative to the wear data of Table 1 are: a^ ¼ 7:101 104 =h;
ni m X X
t i;j expðbti;j Þ t i;j1 expðbti;j1 Þ wi;j expðbti;j Þ expðbt i;j1 Þ i¼1 j¼1 Pm m X t i;n expðbT i Þ ¼0 WðT i Þ Pmi¼1 i i¼1 ½expðbT i Þ 1 i¼1
5. Numerical results
^ i;j ¼ E½w
8 ^ ^ ^ Þ expðbt < c^ ða^ =bÞ½expð bt i;j i;j1 Þ :
0
^0 ^0 b
^0 ^0 b
c^ ½ðt i;j =a Þ ðti;j1 =a Þ
under PLP model under LLP model
From the values of Table 3, the LLP model seems to fit the wear data better than the PLP model. Fig. 4 shows the average of the ‘‘observed’’ wear values at given times obtained by interpolating linearly the measured values, and compares it to the MLE of the expected wear under both the models. The higher accuracy of the LLP model over the PLP model in describing the behavior of the average wear is evident. Fig. 5 shows the measured wear levels of the liners, and depicts the estimated 80% probability interval for W(t), based on the log^ ^ W ¼ ½lnðc^ =E½WðtÞ þ normal approximation with parameters s ^ ^ W ¼ lnðE½WðtÞÞ ^ 2W =2, under the LLP model. The fit s 1Þ1=2 and m of the observed data seems to be very good because the probability band includes a fraction of measured points close to its probability content. The MLE of the parameters of the Weibull thermal crack time distribution of Table 1 are
a^ ¼ 36658 h and b^ ¼ 2:056 (29)
The MLE procedure is preferable to the least-squares procedure based on the Kaplan–Meier estimates (11) as it uses all of the information provided by the survival times. On the contrary, the least-squares procedure does not use the exact value of the survival times but uses, together with the failure times, only the information provided by the number of survival times within two successive ordered failure times.
(31)
The estimated value of b is in agreement with the fatigue mechanism. On the basis of such estimates, the reliability and pdf, Table 3 Comparison of goodness of fit of LLP and PLP models for wearing events.
Estimated log-likelihood Mean squared error
LLP
PLP
134.23 0.0428
136.80 0.0510
ARTICLE IN PRESS
4.5
Measured wear [mm]
4.0 3.5 3.0 2.5 2.0 1.5 Average of observed wear LLP model for wearing events PLP model for wearing events
1.0 0.5 0.0 0
10,000
20,000 30,000 40,000 Operating time [h]
50,000
60,000
Estimated unreliability and failure probabilities
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
1305
1.0 FC (t) 0.8
FC,TC (t)
0.6
0.4
0.2
FC,W (t)
0.0 0
Fig. 4. The average of the observed wear (K) and the MLE of the expected value E[W(t)] under LLP (——) and PLP (- - -) models.
20,000
40,000
60,000
80,000
Operating time [h] Fig. 6. Estimate of liner unreliability FC(t) under the competing risk model (——) and of failure probabilities FC,W(t) and FC,TC(t) (- - -).
5.0 4.0 3.5 3.0 2.5 2.0
1 4 7 10 13 16 19 22 25 28 31
1.5 1.0 0.5 0.0 0
10,000
20,000 30,000 40,000 Operating time [h]
2 5 8 11 14 17 20 23 26 29
50,000
3 6 9 12 15 18 21 24 27 30 33
60,000
Estimated conditional unreliability and failure probabilities
Measured wear [mm]
4.5
Fig. 5. The observed data and the approximate 80% probability interval for W(t) under the LLP model for the wearing events.
1.0 FC (t | t1,W (t1))
0.8 FC,TC (t | t1,W (t1))
0.6 0.4 FC,W (t | t1,W (t1))
0.2 0.0 30,000
35,000
40,000
45,000
50,000
55,000
60,000
Operating time [h] the marginal hazard function, etc., can be easily estimated. For example, the MLE of the expected life, in the absence of the wear mode, is equal to 31,474 h. Using the estimates (30) and (31), all quantities and functions of interest under the competing risk model can be easily estimated. For example, the liner reliability (13) at age t ¼ 40,000 h is equal to 0.255 and the mean lifetime (18) is equal to 29,810 h. Fig. 6 depicts the estimate of the liner unreliability FC(t) under the competing risk model together with the failure probabilities FC,W(t) and FC,TC(t). It shows that the thermal crack is the dominant failure mode and that during the first 20,000 h the wear failure model shows a very low probability of occurrence. The asymptotic values of the failure probabilities FC,W(t) and FC,TC(t) result in 0.129 and 0.871, respectively for wear degradation and thermal cracking, meaning that the wear mechanism is liable for approximately 13% of liner failures. Furthermore, assuming that a liner has not failed by t1 ¼ 30,000 h and the measured wear is W(t1) ¼ 3.5 mm, the conditional reliability (19) at t ¼ 40,000 h, under the competing risk model, is equal to 0.242 and the mean residual lifetime is 7286 h. Fig. 7 shows the estimate of the conditional unreliability F C ðtjt 1 ; Wðt 1 ÞÞ under the competing risk model and the estimate of the conditional failure probabilities F C;W ðtjt 1 ; Wðt 1 ÞÞ and F C;TC ðtjt 1 ; Wðt 1 ÞÞ. It shows that, for such a liner, the thermal crack is the dominant failure mode during the early future 5000 h, where after the wear mechanism becomes dominant. The asymptotic estimates of the conditional failure probabilities of
Fig. 7. Estimate of conditional unreliability under the competing risk model (——) and of conditional failure probabilities (- - -), given t1 ¼ 30,000 h and W(t1) ¼ 3.5 mm.
each failure mode are lim F C;W ðtjt 1 ; Wðt 1 ÞÞ ¼ 0:605 and lim F C;TC ðtjt 1 ; Wðt 1 ÞÞ ¼ 0:395
t!1
t!1
(32) thus showing that, given the survival age t1 ¼ 30,000 h and the accumulated wear W(t1) ¼ 3.5 mm, the failure will more likely be caused by the wear mechanism. However, as depicted in Fig. 8, the marginal hazard function hW(t) for the wear mechanism under the log-linear model (1) shows a non-monotonic behavior: hW(t) increases initially with the liner age up to 60,000 h, and then it decreases. This somewhat anomalous behavior, which mathematically depends on the bounded behavior of the expected number of wearing events of Eq. (2), can find theoretical justification in the so-called strengthening effect that often occurs in degradation processes [10]. Unfortunately, the actual range of the observation period (a few less than 60,000 h) has not allowed the decreasing hW(t) behavior to be possibly observed, and thus the actual presence of a strengthening effect cannot be verified on the basis of the observed data. It is hoped that new wear measurements at greater ages of the liners will be soon available, so that the nonmonotonic behavior of hW(t) could be statistically assessed.
ARTICLE IN PRESS
Marginal hazard functions for wear mode
1306
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
^ By substituting (A.2) into (A.1), the logbut depends on b. likelihood becomes a function of b and c, alone 8 39 2 m P > > > > WðT Þ=c > i < m = 7> 6 X WðT i Þ 6 i¼1 7 ‘ðb; cÞ ¼ ln6 m 7 > c 5> 4P i¼1 > > > ½expðbT i Þ 1 > : ;
1E-4 LLP model for wearing events PLP model for wearing events
9E-5 8E-5 7E-5 6E-5
i¼1
5E-5 þ
4E-5
ni m X X wi;j i¼1 j¼1
3E-5
2E-5
ni m X X
c
ln½expðbt i;j Þ expðbt i;j1 Þ
ln½ðwi;j =cÞ!
i¼1 j¼1
1E-5 0 0
20,000
40,000 60,000 Operating time [h]
80,000
Fig. 8. Comparison of the marginal hazard functions for wear mode under the LLP model (——) and the PLP model (- - -) for the wearing events.
In the light of this uncertainty, an analyst who wishes to be more cautious can use the PLP model to describe the wear events (see, Ref. [6]) and will thus obtain more conservative results. In fact, as shown in Fig. 8, the marginal hazard function hW(t) under the PLP model is a monotonically increasing function with the age t.
6. Conclusions The proposed competing risk model has been shown to give a good fit to the real data of wear degradation and thermal cracking of a number of cylinder liners of Diesel engines, which equip twin ships of the Grimaldi Group, and the obtained results are in agreement with the judgments of Grimaldi experts. The proposed model has not yet been used as input in practical decisions by the Grimaldi management, as the maintenance actions are currently performed by the engine supplier. However, the Grimaldi management intends to develop an alternative maintenance policy that will involve the ships crew. For this reason, they are interested in developing models, such as proposed in this paper, that can enable them to optimize the condition based maintenance activity.
m X WðT i Þ i¼1
c
The first derivative of (A.4) with respect to b is 8
ni m X t i;j exp bt i;j t i;j1 exp bti;j1 @‘ðb; cÞ 1
wi;j ¼ @b c : i¼1 j¼1 exp bt i;j exp bt i;j1 9 m P > > T exp ð bT Þ > i i m = X i¼1 W ðT i Þ m > P i¼1 > expðbT i Þ 1 > ; i¼1
and can be computed without having estimated a and c separately. In particular, from (A.6), the MLE of the sum of the mean wear accumulated during the whole observation period (0, Ti) by all the m liners equals the sum of observed wear levels W(Ti): Pm m m m X X X i¼1 WðT i Þ ^ Þ 1 ¼ ^ ½expðbT WðT i Þ E½WðT i Þ ¼ Pm i ^ i¼1 i¼1 i¼1 ½expðbT i Þ 1 i¼1 (A.7) Considering the random variables wi,j, the LLP that models the wear events is an independent increments process; hence, the wear values wi,j accumulated during non-overlapping intervals are independent random variables with mean and variance, respectively, equal to E½wi;j ¼ ca=b½expðbt i;j Þ expðbt i;j1 Þ and Var½wi;j ¼ cE½wi;j P Thus, the N ¼ m i¼1 ni random variables yk
The log-likelihood function relative to the observed wear data is given by
yk ¼
ni m X X wi;j i¼1 j¼1
c
ni m X X i¼1 j¼1
m aX ½expðbT i Þ 1 b i¼1
wi;j E½wi;j pffiffiffiffiffiffiffiffiffiffiffiffiffi ; E½wi;j
k¼jþ
i1 X
nl
(A.8)
l¼0
have a common mean equal to 0, and common variance equal to c. If the MLE
flnða=bÞ þ ln½expðbti;j Þ expðbt i;j1 Þg
ln½ðwi;j =cÞ!
(A.5)
so that the MLE b^ can be easily obtained by finding the b value that makes the derivative (A.5) null. Such an estimate does not depend on c. From Eqs. (5) and (A.3), the MLE of expected wear accumulated at t is given by Pm WðT i Þ ^ 1 ^ ½expðbtÞ (A.6) E½WðtÞ ¼ Pm i¼1 ^ Þ 1 ½expð bT i i¼1
Appendix A
‘ða; b; cÞ ¼
(A.4)
(A.1)
and by equating to 0 the first derivative of (A.1) with respect to a, we obtain Pm b^ WðT i Þ a^ ¼ Pm i¼1 (A.2) ^ c^ i¼1 ½expðbT i Þ 1 ^ and c^ denote the estimate of a, b, and c, respectively. where a^ , b, From Eq. (A.2), the MLE of y ¼ ca/b is in closed form Pm WðT i Þ (A.3) y^ ¼ c^ a^ =b^ ¼ Pm i¼1 ^ Þ 1 ½expð bT i i¼1
^ Þ expðbt ^ ^ i;j ¼ y^ ½expðbt E½w i;j i;j1 Þ i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni
(A.9)
are used in Eq. (A.8) instead of the true value of E[wi,j], and these estimates are sufficiently accurate, the sample variance of y^ k ¼
^ i;j wi;j E½w qffiffiffiffiffiffiffiffiffiffiffiffiffi ^ i;j E½w
(A.10)
^ i;j are, is close to the elementary wear c. The more accurate E½w and the larger the sample size N, more accurate the estimate of c becomes, provided by the sample variance of y^ k . If the sample mean of y^ k , say Ave½y^ k , is close to 0, then there is confidence that
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
the estimate of c is accurate. Once c and b have been estimated, the estimate of the scale parameter a can be easily obtained from Eq. (A.2).
[7] [8]
References [1] Kodali P, How P, McNulty WD. Methods of improving cylinder liner wear. In: Proceedings of the SAE 2000 world congress, Detroit, 2000. SAE paper 200001-0926. [2] Li S, Csontos AA, Gable BM, Passut CA, Jao T-C. Wear in Cummins M-11/EGR test engine. In: Proceedings of the SAE international fuels & lubricants meeting & exhibition, Reno, 2002. SAE paper 2002-01-1672. [3] Gadow R, Lopez D, Perez J. Carbon deposition reducing coatings for highly loaded large diesel engines. In: Proceedings of the 2003 SAE world congress, Detroit, 2003. SAE paper 2003-01-1100. [4] Sulzer. Service instructions for diesel engine G.M.-Sulzer RT58-Type: 8 RTA 58. Sulzer; 1986. [5] Lu CJ, Meeker WQ. Using degradation measures to estimate a time-to-failure distribution. Technometrics 1993;35(2):161–74. [6] Bocchetti D, Giorgio M, Guida M, Pulcini G. A wear model for condition-based maintenance of cylinder liners in a naval diesel engine. In: Proceedings of the
[9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
1307
8th conference on quality management for organisational and regional development (QMOD), Palermo, 2005. p. 815–25. Crow LH. Reliability analysis for complex systems. In: Proschan F, Serfling RJ, editors. Reliability and biometry. Philadelphia: SIAM; 1974. p. 379–410. Cox DR, Lewis PAW. The statistical analysis of series of events. London: Chapman & Hall; 1966. David HA, Moeschberger ML. The theory of competing risks. London: Charles Griffin & Company; 1978. Gertsbakh IB, Kordonskiy KB. Models of failures. Berlin: Springer; 1969. Weibull W. Fatigue testing and analysis of results. New York: Pergamon Press; 1961. Little RE, Ekvall JC. Statistical analysis of fatigue data. Philadelphia: ASTM; 1979. Murthy DNP, Xie M, Jiang R. Weibull models. New York: Wiley; 2004. Nelson W. Applied life data analysis. New York: Wiley; 1982. Lawless JF. Statistical models and methods for lifetime data. New York: Wiley; 1982. Meeker WQ, Escobar LA. Statistical methods for reliability data. New York: Wiley; 1998. Elandt-Johnson RC, Johnson NL. Survival models and data analysis. New York: Wiley; 1982. David HA, Moeschberger ML. The theory of competing risks. London: Griffin & Company; 1978.
Contents lists available at ScienceDirect
Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress
A competing risk model for the reliability of cylinder liners in marine Diesel engines D. Bocchetti a, M. Giorgio b, M. Guida c, G. Pulcini d, a
Grimaldi Group, Naples, Italy Department of Aerospace and Mechanical Engineering, Second University of Naples, Aversa (NA), Italy Department of Information Engineering and Electrical Engineering, University of Salerno, Fisciano (SA), Italy d Istituto Motori, National Research Council—CNR, Naples, Italy b c
a r t i c l e in fo
abstract
Article history: Received 13 November 2007 Received in revised form 29 January 2009 Accepted 29 January 2009 Available online 10 February 2009
In this paper, a competing risk model is proposed to describe the reliability of the cylinder liners of a marine Diesel engine. Cylinder liners presents two dominant failure modes: wear degradation and thermal cracking. The wear process is described through a stochastic process, whereas the failure time due to the thermal cracking is described by the Weibull distribution. The use of the proposed model allows performing goodness-of-fit test and parameters estimation on the basis of both wear and failure data. Moreover, it enables reliability estimates of the state of the liners to be obtained and the hierarchy of the failure mechanisms to be determined for any given age and wear level of the liner. The model has been applied to a real data set: 33 cylinder liners of Sulzer RTA 58 engines, which equip twin ships of the Grimaldi Group. Estimates of the liner reliability and of other quantities of interest under the competing risk model are obtained, as well as the conditional failure probability and mean residual lifetime, given the survival age and the accumulated wear. Furthermore, the model has been used to estimate the probability that a liner fails due to one of the failure modes when both of these modes act. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Competing risk model Reliability Diesel engine cylinder liners Wear process Thermal cracking
1. Introduction Modern marine Diesel engines must satisfy stringent reliability and availability requirements. Experience shows that the more critical components of heavy-duty Diesel engines are the cylinder liners. These components are very costly and their in-service failures can lead to heavy economic losses. As such, estimating their reliability is crucial when performing long-term economic evaluations and to plan maintenance activities. Unfortunately, reliability estimation of cylinder liners can be difficult to perform due to the scarcity of failure data, and large experimental costs prevent the running of ad hoc reliability tests. As a consequence and for economic reasons, the liners undergo extensive maintenance, aimed at minimizing the number of inservice failures. As such, it is important to make use of models and methodologies that can allow good reliability estimates to be performed based on the available data and all the other technological and experimental information usually possessed by the analyst. Cylinder liners present two dominant failure modes: wear degradation and thermal cracking. In high power marine Diesel
Corresponding author. Tel.: +39 0817177113; fax: +39 081 2396097.
E-mail address: [email protected] (G. Pulcini). 0951-8320/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2009.01.010
propulsion engines, wear usually occurs in the top region of the inner surface of the cylinder liner, where the maximum mechanical and thermal load occurs. Many studies of Diesel engine wear agree that the dominant wear mechanism in this region is due to the high quantity of abrasive particles on the piston surface, generated by the combustion of heavy fuels and oil degradation (soot) (see Refs. [1–3]). The abrasive wear mechanism can be considered as a three-body contact mode. If the lubricant film thickness is less than the soot particle size, then the soot can cause an abrasive action with the liner metal surfaces. Indeed, soot particles are harder than the corresponding engine parts, which further increase the process of abrasion. In addition to abrasive wear, corrosive wear (due to sulphuric acid, nitrous/nitric acids, and water) can also be observed. The thermal cracking of a liner is due to a fatigue mechanism originated by repeated thermal shocks. A thermal shock is caused by a rapid change in the temperature of the cooling fluid that laps the external surface of the liner, and is often aggravated by scale and corrosion in the cooling water spaces as a result of inadequate chemical treatment of the cooling water. These temperature changes occur mainly during the vessel manoeuvre when it is difficult to maintain constant the cooling fluid temperature. As a consequence of the repeated thermal shocks, a small crack can arise in the external surface of the liner and, with time, it propagates towards the inside.
ARTICLE IN PRESS 1300
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
In this paper, a competing risk model is proposed to describe the failure process of a cylinder liner under the hypotheses that (a) the liner fails when the first of the competing failure mechanisms reaches a failure state, (b) the failure mechanisms proceed independently of each other, and (c) each failure mode has a known life distribution. Based on wear measurements taken from a large set of liners at different ages, the wear failure mode is treated as a degradation process for which the failure occurs when a fixed warranty limit is crossed. On the contrary, in the absence of any information on the underlying process of thermal shocks leading to fatigue cracking, the thermal crack mode is treated as a catastrophic failure mode, and the time to failure due to the thermal shocks are modeled via a Weibull distribution. The proposed model has been applied to a real data set referring to the cylinder liners of Sulzer RTA 58 engines, which equip twin ships of the Grimaldi Group. Estimates of the liner reliability and failure rate are obtained based on the competing risk model approach. The probability of failure of the liner is caused by one of the failure modes, when both of these modes act. By its very nature, using the proposed model allows the analyst to: (a) justify the choice of the model on the basis of strong experimental and technological arguments; (b) estimate the reliability of liners using both failure and wear data; (c) use information collected during recurrent maintenance inspections to improve the accuracy of reliability estimates and to perform estimates conditioned to the state of the liners; and (d) recognize the failure modes hierarchy for any given age and wear level of the liner, thus providing the analyst with information about the criticality of the failure mechanisms.
Fig. 1. The cross-section of the cylinder of SULZER RTA 58 engine.
tions are survival times. In three cases (i ¼ 12, 25, and 28), the liner was replaced due to a large accumulated wear, and hence the survival time ti coincides with the time ti;ni of the last inspection. The total cumulated operating time of all the 33 liners was of 913,894 h.
3. Model assumptions 2. The case study
3.1. Independence of failure modes
The failure data set consists of m ¼ 33 cylinder liners of 8cylinder SULZER RTA 58 engines [4], equipping twin ships of the Grimaldi Group. These data were collected from the stage each liner was initially installed, between the dates of January 1, 1994 and February 28, 2007. The liners were recurrently inspected and the accumulated wear was measured. Measurements were taken by positioning a caliper inside the cylinder bore at the position of the Top Dead Center (TDC) (Fig. 1): the point of maximum wear [2]. The sensitivity of the caliper was 0.05 mm, and as such each measure was rounded up to the nearest 0.05. Furthermore, when a liner failed due to thermal cracking, it was replaced with a new one. Thus, the wear data set consisted of the number ni (i ¼ 1, y, 33) of inspections performed on each liner, the age ti,j (j ¼ 1, y, ni) (in operating hours) of the liner at the time of inspection j, and the wear Wi,jW(ti,j) (in mm) measured during P each inspection (see Table 1). A total of N ¼ m i¼1 ni ¼ 64 measures of wear were recorded. The liner i ¼ 32 failed due to thermal cracking before the first inspection was performed, and thus no wear data was subsequently recorded from it. Thermal shocks were not recorded, even if the cooling fluid temperature was monitored, such that the data concerning the thermal crack (given in Table 1) consist only of the age ti (i ¼ 1, y, 33) at which either the failure occurred (di ¼ 1) or did not occur (di ¼ 0), due to thermal shock. A total of 22 failures by thermal cracking were observed, and the remaining 11 observa-
The two failures mechanisms, namely the wear and the thermal crack, proceed independently of each other (at least until a failure occurs). This assumption is justified by the fact that the physical mechanisms that cause or influence the evolution of each failure mode are quite different and unrelated. In fact, wear that occurs on the inner surface of the liner cannot affect the mechanical strength of the liner. Indeed, the maximum admissible wear (4 mm) is negligible with respect to the initial liner thickness (100 mm). Furthermore, the point at which maximum wear occurs is far enough from the point at which thermal crack initiation occurs and subsequently propagates towards the inside. Also, the changes in the temperature of the cooling fluid, which are the driving mechanism for thermal cracking, do not influence the wear formation process, the size of abrasive particles, or lubrication oil degradation. The independence assumption implies that the probability of a liner failure in (t,t+dt) from one of the two failure mechanisms, given survival to time t, is the same as to whether that mechanism were the only failure mechanism acting on the liner. 3.2. The wear model Using the wear measurements taken on the liners at different ages, the wear failure mode is treated as a degradation process for
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
1301
Table 1 Wear and thermal crack data of cylinder liners. i
ni
ti,1 (h)
Wi,1 (mm)
ti,2 (h)
Wi,2 (mm)
ti,3 (h)
Wi,3 (mm)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
3 2 2 2 3 3 3 3 3 4 3 4 3 2 1 1 1 1 1 1 1 1 1 1 4 1 1 3 2 1 1 0 1
11,300 11,360 11,300 12,300 14,810 9700 10,000 6860 2040 7540 8510 18,320 10,000 9350 13,200 7700 7700 8250 3900 7700 9096 19,800 10,450 12,100 12,000 8800 2200 33,000 8800 8250 18,755
0.90 0.80 1.50 1.00 1.90 1.10 1.20 0.50 0.40 0.50 0.80 2.20 2.10 0.85 2.00 1.05 1.60 0.90 1.15 1.20 0.50 1.60 0.40 1.00 1.95 1.40 0.40 2.90 0.50 0.70 1.15
14,680 17,200 21,970 16,300 18,700 19,710 30,450 17,200 12,580 8840 14,930 25,310 16,620 15,970
1.30 1.35 2.00 1.35 2.25 2.60 2.75 1.45 2.00 1.10 1.45 3.00 2.75 1.20
31,270
2.85
28,000 30,450 37,310 24,710 16,620 9770 21,560 37,310 30,000
2.75 3.00 3.05 2.15 2.35 1.15 1.90 3.70 3.60
27,300
2.70
49,500
3.15
38,500 27,500
3.25 2.15
55,460
4.10
8490
0.95
which the degradation warranty limit Wmax ¼ 4 mm is considered [5]. The wear process is then modeled through a point stochastic process, by assuming (as suggested in Ref. [6]) that: (a) the wear of a cylinder liner at a given age t is the accumulation of isolated damage (elementary wear) of equal magnitude c that has occurred by age t; (b) the events that cause an elementary wear occur randomly in time; and (c) the probability that damage occurs in a given time interval is independent of the previously accumulated damage.
0
0
a0 40; 0ob o1.
However, since the wear process cannot grow unlimitedly, a stochastic process that bounds the expected number of wearing events offers a more realistic modeling approach. Considering this, the wearing process can now be described through a nonhomogeneous Poisson process (NHPP) with intensity function:
lðtÞ ¼ a expðbtÞ;
a40; 1obo1
(1)
often referred to as log-linear process (LLP) [8]. Based on the LLP assumption, the expected number of wear events up to the operating time t, say MðtÞ E½NðtÞ, is given by MðtÞ E½NðtÞ ¼ ða=bÞ½expðbtÞ 1
(2)
Note that, when bo0, M(t) tends to a/b as the age t of the liner increases. Assumptions are made in terms of an independent increments process and of constant elementary damages in order to avoid
Wi,4 (mm)
16,300
2.10
45,000
3.95
56,120
4.05
ti (h)
di
36,370 28,930 27,970 21,830 39,450 39,500 42,320 25,200 27,750 25,680 29,900 45,000 37,860 25,100 18,270 18,650 26,770 16,870 11,600 14,300 14,596 31,900 25,300 23,720 56,120 20,950 16,738 55,460 28,100 8250 31,330 5430 16,790
1 1 1 1 0 1 0 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 0 0 0 1 0 1 0 1 1 1
that some estimation issues, such as those discussed later in Section 4.1, become unsolvable. However, despite these simplifications, the proposed model proved to be suitable for this case study. On the basis of the above assumptions, the wear accumulated at the operating time t is modeled by the stochastic process: WðtÞ ¼ cNðtÞ
(3)
so that the probability distribution of the wear w(t1, t2) accumulated in the time interval (t1, t2) is P½wðt1 ; t 2 Þ ¼ w ¼ P½Nðt 1 ; t 2 Þ ¼ w=c
Ref. [6] describes the wearing events process N(t) through a power law process (PLP) [7] with decreasing intensity function:
l0 ðtÞ ¼ ðb0 =a0 Þðt=a0 Þb 1 ;
ti,4 (h)
¼
½Mðt 2 Þ Mðt 1 Þw=c exp½Mðt 2 Þ þ Mðt 1 Þ ðw=cÞ!
(4)
where N(t1, t2) denotes the number of wear events in the interval (t1, t2) and w ¼ 0, c, 2c,y . For each NHPP, Var[N(t)] ¼ M(t). Therefore, the mean and variance of the wear accumulated until t are, respectively E½WðtÞ ¼ cMðtÞ and Var½WðtÞ ¼ c2 MðtÞ
(5)
Since the time at which the wear level crosses the warranty limit Wmax ¼ 4 mm is considered as a failure time, the wear mode reliability of a liner, defined as the hypothetical probability that, in absence of the other failure modes, the wear accumulated until t does not exceed the 4 mm limit [9], is given by RW ðtÞ ¼ P½WðtÞpW max ¼ P½NðtÞpK max ¼
K max X
½MðtÞk exp½MðtÞ ¼ 1 IG½MðtÞ; K max þ 1 k! k¼0
(6)
where K max ¼ int½W max =c is the greatest integer equal to or less Ry than Wmax/c, and IGðy; kÞ ¼ 0 xk1 expðxÞdx=GðkÞ is the incom-
ARTICLE IN PRESS 1302
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
plete Gamma function. The probability density function (pdf) is f W ðtÞ ¼
a½MðtÞ
GðK max þ 1Þ
exp½bt MðtÞ;
t40
(7)
RW ðtjWðt 1 ÞÞ ¼ P½WðtÞpW max jWðt 1 Þ ¼ P½Nðt 1 ; tÞpDK max DX K max k¼0
½MðtÞ Mðt 1 Þk exp½MðtÞ þ Mðt 1 Þ k!
¼ 1 IG½MðtÞ Mðt 1 Þ; DK max þ 1
(8)
where DK max ¼ int½fW max Wðt 1 Þg=c is the greatest integer equal to or less than {WmaxW(t1)}/c. Fig. 2 shows the measured wear accumulated by the cylinder liners, where the measured points for each liner have been linearly interpolated for graphical representation. The wear behavior is in keeping with the literature results (see, for example, Ref. [10]) and the technological characteristics: the early phase is a wear accommodation period in which the rate of wear decreases quickly, followed by a phase in which the rate decreases more slowly. In the present application, the accommodation phase occupies the large part of the entire observation time, since the warranty clauses impose that a cylinder liner is replaced before it accumulates a relatively ‘‘low’’ wear level (4 mm), in order to avoid potentially catastrophic and very expensive failures. Moreover, the observed (interpolated) processes are closely interwoven, thus indicating that the rate of wear of an individual liner changes continually, and that these changes are of a random nature. 3.3. The thermal crack model In absence of information regarding the underlying process of thermal shocks, leading to fatigue cracking, the thermal crack mode is treated as a hard (or catastrophic) failure mode, and its reliability at the age t (i.e. the hypothetical probability that a thermal crack failure occurs after t in absence of the wear mode) is assumed to be Weibull RTC ðtÞ ¼ exp½ðt=aÞb ;
a; b40
(9)
This assumption is justified on the basis of theoretical considerations since fatigue life data is often shown to be adequately analyzed using the Weibull distribution (see Refs. [11–13]), and is supported by graphical analysis. In particular, the graphical analysis was performed by plotting, on Weibull paper, the Kaplan–Meier estimate [14] of the reliability
Measured wear [mm]
4.0 3.5 3.0 2.5 1 4 7 10 13 16 19 22 25 28 31
2.0 1.5 1.0 0.5
2 5 8 11 14 17 20 23 26 29
3 6 9 12 15 18 21 24 27 30 33
0.0 0
R¯ i R¯ TC ðti Þ ¼
Y rl 1
10,000
20,000
30,000
40,000
Operating time [h] Fig. 2. The observed wear processes.
50,000
60,000
(10)
rl
l2Si
Given the wear W(t1)oWmax accumulated by time t1, the conditional probability that the warranty limit is not exceeded before t4t1 is given by
¼
at any catastrophic failure time ti (i:di ¼ 1) versus ti:
K max
where Si is the set of liners whose catastrophic failure occurred at tlpti, and rl is the reverse rank, i.e., the number of failure or censoring times no less than rl (see Table 2, where the catastrophic failure and censoring times were ordered for convenience reason). For example, S4 ¼ (4, 15, 16, 18, 19, 20, 21, 27, 32, 33) and r4 ¼ 23. As depicted in the Weibull plot of Fig. 3, the points roughly follow a straight line and give no obvious evidence that the thermal crack data do not fit a Weibull distribution. Furthermore, the high determination coefficient, say R2 ¼ 0.95, shows a high goodness of fit. Table 2 also gives both the estimated standard deviation of the estimated reliability [15], based on the following function:
sfR¯ i g ¼ R¯ i
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ½r l ðr l 1Þ1
(11)
l2Si
and the approximate g ¼ 0.80 confidence interval on RTC(ti), which is based on logit transformation [16] and is given by
R¯ i R¯ i ; R¯ i þ ð1 R¯ i Þs R¯ i þ ð1 R¯ i Þ=s
(12)
where s ¼ expfz1g=2 sfR¯ i g=½R¯ i ð1 R¯ i Þg and z1g/2 is the (1g/2) quantile of the Normal standard distribution. As shown in Fig. 3, the fitted line falls within the 80% confidence intervals. Note that the endpoints of Eq. (12) always lie between 0 and 1. Table 2 Kaplan–Meier estimate and 80% confidence interval of reliability for thermal crack mode. Liner i
ti
di
ri
R¯ i
sfR¯ i g
80% confidence interval
32 19 20 21 27 33 18 15 16 26 4 24 14 8 23 10 17 9 3 30 29 2 11 31 22 1 13 5 6 7 12 28 25
5430 11,600 14,300 14,596 16,738 16,790 16,870 18,270 18,650 20,950 21,830 23,720 25,100 25,200 25,300 25,680 26,770 27,750 27,970 28,050 28,100 28,930 29,990 31,330 31,900 36,370 37,860 39,450 39,500 42,320 45,000 55,460 56,120
1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 1 0 1 1 0 1 1 1 1 1 1 0 0 1 0 0 0 0
33 32 31 30 29 28 27 26 25
0.9697 0.9394 0.9091 0.8788 0.8485 0.8182 0.7879 0.7576 0.7273
0.0298 0.0298 0.0415 0.0500 0.0567 0.0622 0.0669 0.0709 0.0742
(0.8970, 0.9916) (0.8879, 0.9681) (0.8401, 0.9501) (0.7990, 0.9297) (0.7609, 0.9079) (0.7247, 0.8850) (0.6898, 0.8612) (0.6559, 0.8367) (0.6228, 0.8116)
23
0.6957
0.0774
(0.5885, 0.7851)
20 19 18
0.6609 0.6261 0.5913
0.0810 0.0839 0.0861
(0.5508, 0.7559) (0.5141, 0.7260) (0.4782, 0.6955)
16 15
0.5543 0.5174
0.0883 0.0898
(0.4403, 0.6629) (0.4034, 0.6296)
13 12 11 10 9 8
0.4776 0.4378 0.3980 0.3582 0.3184 0.2786
0.0913 0.0920 0.0918 0.0909 0.0890 0.0864
(0.3638, 0.5937) (0.3254, 0.5569) (0.2880, 0.5193) (0.2517, 0.4808) (0.2164, 0.4414) (0.1821, 0.4011)
5
0.2229
0.0852
(0.1325, 0.3501)
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
under the competing risk model, is
0.80 0.70 0.60 0.50 0.40
hC ðtÞ ¼ hW ðtÞ þ hTC ðtÞ ¼ b1
þ
0.30
a½MðtÞkmax exp½bt MðtÞ Gðkmax þ 1Þf1 IG½MðtÞ; kmax þ 1g
b t
(17)
a a
The mean lifetime can be obtained by numerical integration of the liner reliability (13): Z 1 RC ðtÞdt (18) EC ½T ¼
0.20 1 - RTC (t)
1303
0
0.10
Given that the thermal crack has not occurred until t1 and the accumulated wear at t1 is Wðt 1 ÞoW max , the conditional probability that the liner does not fail before t4t1 is given by " b # t1 b t RC ðtjt 1 ; Wðt 1 ÞÞ ¼ f1 IG½MðtÞ; Dkmax þ 1gexp
0.05
a
0.02
00 ,0
00 ,0
The mean residual life, given the survival age t1 and the accumulated wear W(t1), is given by the integration of (19) over the interval (t1,N). Furthermore, the conditional pdfs at tXt1 for the wear and thermal crack mode are, respectively
50
Operating time (in log scale)
f C;W ðtjt 1 ; Wðt 1 ÞÞ ¼
Fig. 3. Graphical analysis of thermal crack data on the Weibull plot.
a½MðtÞ Mðt 1 ÞDkmax GðDkmax þ 1Þ " exp bt MðtÞ þ Mðt 1 Þ
3.4. The competing risk model
f C;TC ðtjt 1 ; Wðt 1 ÞÞ ¼
Under the assumption of independence, the reliability, under the competing risk model, is the product of the failure mode reliability functions (6) and (9) [9,17,18]: RC ðtÞ ¼ RW ðtÞRTC ðtÞ
(13)
The probability density functions of each failure mode in the presence of both the modes are given by f C;W ðtÞ ¼ f W ðtÞRTC ðtÞ ¼
a½MðtÞkmax
Gðkmax þ 1Þ b t b1
exp½bt MðtÞ ðt=aÞb
f C;TC ðtÞ ¼ f TC ðtÞRW ðtÞ ¼ a a " # b t f1 IG½MðtÞ; kmax þ 1g exp
a
f C ðtÞ ¼ f C;W ðtÞ þ f C;TC ðtÞ
(14)
(15)
Thus, from a computational viewpoint, it is more convenient to obtain first the failure probability FC,W(t) by integrating numerically its pdf fC,W(t), and then derive FC,TC(t) from F C;TC ðtÞ ¼ 1 RW ðtÞRTC ðtÞ F C;W ðtÞ
and b t b1
a a
b t
a
þ
b # t1
a
" b # t1 b t exp
a
a
f1 IG½MðtÞ Mðt 1 Þ; Dkmax þ 1g
(20)
The conditional failure probabilities F C;W ðtjt 1 ; Wðt 1 ÞÞ and F C;TC ðtjt 1 ; Wðt 1 ÞÞ are then given by the (numerical) integration of (20). Their asymptotic value, for t-N, provides the probability that wear degradation and thermal cracking failure occur, respectively, given the state of the liner at time t1.
4. Estimation procedure
The associated failure probabilities, say FC,W(t) and FC,TC(t), are given by the (numerical) integration of (14), and their asymptotic value, for t-N, is less than 1 and provides the proportion of failures caused by each failure mode during the whole life of a liner. The liner pdf and cumulative distribution function (Cdf), under the competing risk model, are thus given by the following, respectively:
F C ðtÞ ¼ F C;W ðtÞ þ F C;TC ðtÞ
a
(19)
40
30
,0
00 ,0 20
10
5,
,0
00
00
0
0.01
00
Approx. 80% Confidence Interval Non Parametric Estimate Weibull Fit
(16)
Under the independence assumption, the failure mode hazard functions hC,W(t) and hC,TC(t) equal the marginal hazard functions hW(t) and hTC(t), respectively, so that the liner hazard function,
Under the assumption of Section 3, the likelihood function relative to the observed data is 9 8 ni m Y
i¼1 m b X ti i¼1
#
a
i¼1
a a
(21)
where wi;j ¼ Wðt i;j Þ Wðt i;j1 Þ is the wear accumulated by the liner i during the inspection interval (ti,j1, ti,j) (when W(ti,0)0 and ti,00), and T i t i;ni is the age of the liner i at the time of the last inspection. Since the stochastic models for the two failure modes are indexed by uncommon parameters, then the parameter estimation can be performed separately for each mode. 4.1. Estimate of wear model parameters For the wear model (3), both the parameters a and c of the intensity function (1) act as multiplicative factors for the accumulated wear. Hence, the maximum likelihood estimation (MLE) procedure cannot be used to estimate a and c individually.
ARTICLE IN PRESS 1304
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
Thus, to overcome this estimation problem, the following fivestep estimation procedure, that combines MLE and first moments estimation methods, is used (see Appendix A for the mathematical proof): 1. Obtain the MLE b^ by solving numerically the equation
(22)
3. Estimate the expected wear accumulated by each liner i during each interval (ti,j1, ti,j) i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni (24) and compute the N ¼ y^ k ¼
^ i;j wi;j E½w qffiffiffiffiffiffiffiffiffiffiffiffiffi ; ^ i;j E½w
Pm
i¼1 ni
values of
i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni ; k ¼ j þ
i1 X
nl
l¼0
(25) 4. Estimate the elementary wear c as PN ðy^ Ave½y^ k Þ2 c^ ¼ Var½y^ k ¼ k¼1 k N1 P ^ ^ where Ave½y^ k ¼ N k¼1 yk =N is the sample mean of yk . 5. Estimate the LLP parameter a as ^ c^ a^ ¼ y^ b=
(26)
(27)
Note that, since E½WðtÞ ¼ y½expðbtÞ 1, its estimate is a MLE. 4.2. Estimate of thermal crack model parameters For the thermal crack model (9), the MLE procedure for Weibull randomly censored sample [15] is used. In particular, the MLE of the parameter b is the numerical solution of Pm b Pm ti lnðti Þ 1 d lnðti Þ þ i¼1 i¼1 ¼0 (28) Pmi Pm b i¼1 di
t
i¼1 i
and the MLE of a is 2P
a^ ¼
31=b^ b^ m i 5 4 Pi¼1 m i¼1 di
t
c^ ¼ 0:1810 mm
The MLE of b is negative, thus indicating that the rate of wear decreases with the liner age. The sample mean of y^ k is equal to 0.0351. This value is much less than the sample standard deviation sðy^ k Þ ¼ 0:4254, thus giving confidence that the estimate of c is sufficiently accurate. ^ and c^ have been obtained without The estimates a^ , b, constraining the ratios wi,j/c, which give the number of wear events experienced by the liner i during the inspection interval (ti,j1, ti,j), to be integer values. Such a lack of restriction does not imply any mathematical issues, even though Eq. (A.1) no longer represents a probability statement. However, because the isolated damages are not of equal magnitude and the data set consists of rounded measures of wear, the lack of restriction on the ratios wi,j/ c does not produce misleading estimates. For a comparison, the wear data are analyzed also within the 0 0 0 PLP model with intensity function l ðtÞ ¼ ðb =a0 Þðt=a0 Þb 1 proposed in Ref. [6]. The resulting MLE of the PLP parameters are 0 0 a^ ¼ 776:5 h and b^ ¼ 0:7458, and the estimate of the elementary 0 wear is c^ ¼ 0:1735 mm. The goodness of fit of the wear models is compared through: (a) the estimated log-likelihoods relative to ^ c^ Þ and lða^ 0 ; b^ 0 ; c^ 0 Þ, and (b) the mean squared the wear data, say lða^ ; b; error (MSE), defined as Pm Pni 2 ^ i¼1 j¼1 ðwi;j E½wi;j Þ MSE ¼ m where
The accuracy of the estimate of c, and hence of a, depends ^ i;j and on the sample size N. A larger N on the accuracy of E½w provides a more accurate estimate of c, as given by Eq. (26). If the sample mean of y^ k is close to 0, then we can be confident that the estimate of c is accurate.
b
b^ ¼ 2:423 105 =h;
(30)
where W(Ti) is the accumulated wear of the liner i by the last inspection. 2. Obtain the MLE of y ¼ ca/b from Pm WðT i Þ (23) y^ ¼ Pm i¼1 ^ i¼1 ½expðbT i Þ 1
^ Þ expðbt ^ ^ i;j ¼ y^ ½expðbt E½w i;j i;j1 Þ;
On the basis of the estimation procedure of Section 4.1, the estimate of the wear process parameters relative to the wear data of Table 1 are: a^ ¼ 7:101 104 =h;
ni m X X
t i;j expðbti;j Þ t i;j1 expðbti;j1 Þ wi;j expðbti;j Þ expðbt i;j1 Þ i¼1 j¼1 Pm m X t i;n expðbT i Þ ¼0 WðT i Þ Pmi¼1 i i¼1 ½expðbT i Þ 1 i¼1
5. Numerical results
^ i;j ¼ E½w
8 ^ ^ ^ Þ expðbt < c^ ða^ =bÞ½expð bt i;j i;j1 Þ :
0
^0 ^0 b
^0 ^0 b
c^ ½ðt i;j =a Þ ðti;j1 =a Þ
under PLP model under LLP model
From the values of Table 3, the LLP model seems to fit the wear data better than the PLP model. Fig. 4 shows the average of the ‘‘observed’’ wear values at given times obtained by interpolating linearly the measured values, and compares it to the MLE of the expected wear under both the models. The higher accuracy of the LLP model over the PLP model in describing the behavior of the average wear is evident. Fig. 5 shows the measured wear levels of the liners, and depicts the estimated 80% probability interval for W(t), based on the log^ ^ W ¼ ½lnðc^ =E½WðtÞ þ normal approximation with parameters s ^ ^ W ¼ lnðE½WðtÞÞ ^ 2W =2, under the LLP model. The fit s 1Þ1=2 and m of the observed data seems to be very good because the probability band includes a fraction of measured points close to its probability content. The MLE of the parameters of the Weibull thermal crack time distribution of Table 1 are
a^ ¼ 36658 h and b^ ¼ 2:056 (29)
The MLE procedure is preferable to the least-squares procedure based on the Kaplan–Meier estimates (11) as it uses all of the information provided by the survival times. On the contrary, the least-squares procedure does not use the exact value of the survival times but uses, together with the failure times, only the information provided by the number of survival times within two successive ordered failure times.
(31)
The estimated value of b is in agreement with the fatigue mechanism. On the basis of such estimates, the reliability and pdf, Table 3 Comparison of goodness of fit of LLP and PLP models for wearing events.
Estimated log-likelihood Mean squared error
LLP
PLP
134.23 0.0428
136.80 0.0510
ARTICLE IN PRESS
4.5
Measured wear [mm]
4.0 3.5 3.0 2.5 2.0 1.5 Average of observed wear LLP model for wearing events PLP model for wearing events
1.0 0.5 0.0 0
10,000
20,000 30,000 40,000 Operating time [h]
50,000
60,000
Estimated unreliability and failure probabilities
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
1305
1.0 FC (t) 0.8
FC,TC (t)
0.6
0.4
0.2
FC,W (t)
0.0 0
Fig. 4. The average of the observed wear (K) and the MLE of the expected value E[W(t)] under LLP (——) and PLP (- - -) models.
20,000
40,000
60,000
80,000
Operating time [h] Fig. 6. Estimate of liner unreliability FC(t) under the competing risk model (——) and of failure probabilities FC,W(t) and FC,TC(t) (- - -).
5.0 4.0 3.5 3.0 2.5 2.0
1 4 7 10 13 16 19 22 25 28 31
1.5 1.0 0.5 0.0 0
10,000
20,000 30,000 40,000 Operating time [h]
2 5 8 11 14 17 20 23 26 29
50,000
3 6 9 12 15 18 21 24 27 30 33
60,000
Estimated conditional unreliability and failure probabilities
Measured wear [mm]
4.5
Fig. 5. The observed data and the approximate 80% probability interval for W(t) under the LLP model for the wearing events.
1.0 FC (t | t1,W (t1))
0.8 FC,TC (t | t1,W (t1))
0.6 0.4 FC,W (t | t1,W (t1))
0.2 0.0 30,000
35,000
40,000
45,000
50,000
55,000
60,000
Operating time [h] the marginal hazard function, etc., can be easily estimated. For example, the MLE of the expected life, in the absence of the wear mode, is equal to 31,474 h. Using the estimates (30) and (31), all quantities and functions of interest under the competing risk model can be easily estimated. For example, the liner reliability (13) at age t ¼ 40,000 h is equal to 0.255 and the mean lifetime (18) is equal to 29,810 h. Fig. 6 depicts the estimate of the liner unreliability FC(t) under the competing risk model together with the failure probabilities FC,W(t) and FC,TC(t). It shows that the thermal crack is the dominant failure mode and that during the first 20,000 h the wear failure model shows a very low probability of occurrence. The asymptotic values of the failure probabilities FC,W(t) and FC,TC(t) result in 0.129 and 0.871, respectively for wear degradation and thermal cracking, meaning that the wear mechanism is liable for approximately 13% of liner failures. Furthermore, assuming that a liner has not failed by t1 ¼ 30,000 h and the measured wear is W(t1) ¼ 3.5 mm, the conditional reliability (19) at t ¼ 40,000 h, under the competing risk model, is equal to 0.242 and the mean residual lifetime is 7286 h. Fig. 7 shows the estimate of the conditional unreliability F C ðtjt 1 ; Wðt 1 ÞÞ under the competing risk model and the estimate of the conditional failure probabilities F C;W ðtjt 1 ; Wðt 1 ÞÞ and F C;TC ðtjt 1 ; Wðt 1 ÞÞ. It shows that, for such a liner, the thermal crack is the dominant failure mode during the early future 5000 h, where after the wear mechanism becomes dominant. The asymptotic estimates of the conditional failure probabilities of
Fig. 7. Estimate of conditional unreliability under the competing risk model (——) and of conditional failure probabilities (- - -), given t1 ¼ 30,000 h and W(t1) ¼ 3.5 mm.
each failure mode are lim F C;W ðtjt 1 ; Wðt 1 ÞÞ ¼ 0:605 and lim F C;TC ðtjt 1 ; Wðt 1 ÞÞ ¼ 0:395
t!1
t!1
(32) thus showing that, given the survival age t1 ¼ 30,000 h and the accumulated wear W(t1) ¼ 3.5 mm, the failure will more likely be caused by the wear mechanism. However, as depicted in Fig. 8, the marginal hazard function hW(t) for the wear mechanism under the log-linear model (1) shows a non-monotonic behavior: hW(t) increases initially with the liner age up to 60,000 h, and then it decreases. This somewhat anomalous behavior, which mathematically depends on the bounded behavior of the expected number of wearing events of Eq. (2), can find theoretical justification in the so-called strengthening effect that often occurs in degradation processes [10]. Unfortunately, the actual range of the observation period (a few less than 60,000 h) has not allowed the decreasing hW(t) behavior to be possibly observed, and thus the actual presence of a strengthening effect cannot be verified on the basis of the observed data. It is hoped that new wear measurements at greater ages of the liners will be soon available, so that the nonmonotonic behavior of hW(t) could be statistically assessed.
ARTICLE IN PRESS
Marginal hazard functions for wear mode
1306
D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
^ By substituting (A.2) into (A.1), the logbut depends on b. likelihood becomes a function of b and c, alone 8 39 2 m P > > > > WðT Þ=c > i < m = 7> 6 X WðT i Þ 6 i¼1 7 ‘ðb; cÞ ¼ ln6 m 7 > c 5> 4P i¼1 > > > ½expðbT i Þ 1 > : ;
1E-4 LLP model for wearing events PLP model for wearing events
9E-5 8E-5 7E-5 6E-5
i¼1
5E-5 þ
4E-5
ni m X X wi;j i¼1 j¼1
3E-5
2E-5
ni m X X
c
ln½expðbt i;j Þ expðbt i;j1 Þ
ln½ðwi;j =cÞ!
i¼1 j¼1
1E-5 0 0
20,000
40,000 60,000 Operating time [h]
80,000
Fig. 8. Comparison of the marginal hazard functions for wear mode under the LLP model (——) and the PLP model (- - -) for the wearing events.
In the light of this uncertainty, an analyst who wishes to be more cautious can use the PLP model to describe the wear events (see, Ref. [6]) and will thus obtain more conservative results. In fact, as shown in Fig. 8, the marginal hazard function hW(t) under the PLP model is a monotonically increasing function with the age t.
6. Conclusions The proposed competing risk model has been shown to give a good fit to the real data of wear degradation and thermal cracking of a number of cylinder liners of Diesel engines, which equip twin ships of the Grimaldi Group, and the obtained results are in agreement with the judgments of Grimaldi experts. The proposed model has not yet been used as input in practical decisions by the Grimaldi management, as the maintenance actions are currently performed by the engine supplier. However, the Grimaldi management intends to develop an alternative maintenance policy that will involve the ships crew. For this reason, they are interested in developing models, such as proposed in this paper, that can enable them to optimize the condition based maintenance activity.
m X WðT i Þ i¼1
c
The first derivative of (A.4) with respect to b is 8
ni m X t i;j exp bt i;j t i;j1 exp bti;j1 @‘ðb; cÞ 1
wi;j ¼ @b c : i¼1 j¼1 exp bt i;j exp bt i;j1 9 m P > > T exp ð bT Þ > i i m = X i¼1 W ðT i Þ m > P i¼1 > expðbT i Þ 1 > ; i¼1
and can be computed without having estimated a and c separately. In particular, from (A.6), the MLE of the sum of the mean wear accumulated during the whole observation period (0, Ti) by all the m liners equals the sum of observed wear levels W(Ti): Pm m m m X X X i¼1 WðT i Þ ^ Þ 1 ¼ ^ ½expðbT WðT i Þ E½WðT i Þ ¼ Pm i ^ i¼1 i¼1 i¼1 ½expðbT i Þ 1 i¼1 (A.7) Considering the random variables wi,j, the LLP that models the wear events is an independent increments process; hence, the wear values wi,j accumulated during non-overlapping intervals are independent random variables with mean and variance, respectively, equal to E½wi;j ¼ ca=b½expðbt i;j Þ expðbt i;j1 Þ and Var½wi;j ¼ cE½wi;j P Thus, the N ¼ m i¼1 ni random variables yk
The log-likelihood function relative to the observed wear data is given by
yk ¼
ni m X X wi;j i¼1 j¼1
c
ni m X X i¼1 j¼1
m aX ½expðbT i Þ 1 b i¼1
wi;j E½wi;j pffiffiffiffiffiffiffiffiffiffiffiffiffi ; E½wi;j
k¼jþ
i1 X
nl
(A.8)
l¼0
have a common mean equal to 0, and common variance equal to c. If the MLE
flnða=bÞ þ ln½expðbti;j Þ expðbt i;j1 Þg
ln½ðwi;j =cÞ!
(A.5)
so that the MLE b^ can be easily obtained by finding the b value that makes the derivative (A.5) null. Such an estimate does not depend on c. From Eqs. (5) and (A.3), the MLE of expected wear accumulated at t is given by Pm WðT i Þ ^ 1 ^ ½expðbtÞ (A.6) E½WðtÞ ¼ Pm i¼1 ^ Þ 1 ½expð bT i i¼1
Appendix A
‘ða; b; cÞ ¼
(A.4)
(A.1)
and by equating to 0 the first derivative of (A.1) with respect to a, we obtain Pm b^ WðT i Þ a^ ¼ Pm i¼1 (A.2) ^ c^ i¼1 ½expðbT i Þ 1 ^ and c^ denote the estimate of a, b, and c, respectively. where a^ , b, From Eq. (A.2), the MLE of y ¼ ca/b is in closed form Pm WðT i Þ (A.3) y^ ¼ c^ a^ =b^ ¼ Pm i¼1 ^ Þ 1 ½expð bT i i¼1
^ Þ expðbt ^ ^ i;j ¼ y^ ½expðbt E½w i;j i;j1 Þ i ¼ 1; . . . ; m; j ¼ 1; . . . ; ni
(A.9)
are used in Eq. (A.8) instead of the true value of E[wi,j], and these estimates are sufficiently accurate, the sample variance of y^ k ¼
^ i;j wi;j E½w qffiffiffiffiffiffiffiffiffiffiffiffiffi ^ i;j E½w
(A.10)
^ i;j are, is close to the elementary wear c. The more accurate E½w and the larger the sample size N, more accurate the estimate of c becomes, provided by the sample variance of y^ k . If the sample mean of y^ k , say Ave½y^ k , is close to 0, then there is confidence that
ARTICLE IN PRESS D. Bocchetti et al. / Reliability Engineering and System Safety 94 (2009) 1299–1307
the estimate of c is accurate. Once c and b have been estimated, the estimate of the scale parameter a can be easily obtained from Eq. (A.2).
[7] [8]
References [1] Kodali P, How P, McNulty WD. Methods of improving cylinder liner wear. In: Proceedings of the SAE 2000 world congress, Detroit, 2000. SAE paper 200001-0926. [2] Li S, Csontos AA, Gable BM, Passut CA, Jao T-C. Wear in Cummins M-11/EGR test engine. In: Proceedings of the SAE international fuels & lubricants meeting & exhibition, Reno, 2002. SAE paper 2002-01-1672. [3] Gadow R, Lopez D, Perez J. Carbon deposition reducing coatings for highly loaded large diesel engines. In: Proceedings of the 2003 SAE world congress, Detroit, 2003. SAE paper 2003-01-1100. [4] Sulzer. Service instructions for diesel engine G.M.-Sulzer RT58-Type: 8 RTA 58. Sulzer; 1986. [5] Lu CJ, Meeker WQ. Using degradation measures to estimate a time-to-failure distribution. Technometrics 1993;35(2):161–74. [6] Bocchetti D, Giorgio M, Guida M, Pulcini G. A wear model for condition-based maintenance of cylinder liners in a naval diesel engine. In: Proceedings of the
[9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
1307
8th conference on quality management for organisational and regional development (QMOD), Palermo, 2005. p. 815–25. Crow LH. Reliability analysis for complex systems. In: Proschan F, Serfling RJ, editors. Reliability and biometry. Philadelphia: SIAM; 1974. p. 379–410. Cox DR, Lewis PAW. The statistical analysis of series of events. London: Chapman & Hall; 1966. David HA, Moeschberger ML. The theory of competing risks. London: Charles Griffin & Company; 1978. Gertsbakh IB, Kordonskiy KB. Models of failures. Berlin: Springer; 1969. Weibull W. Fatigue testing and analysis of results. New York: Pergamon Press; 1961. Little RE, Ekvall JC. Statistical analysis of fatigue data. Philadelphia: ASTM; 1979. Murthy DNP, Xie M, Jiang R. Weibull models. New York: Wiley; 2004. Nelson W. Applied life data analysis. New York: Wiley; 1982. Lawless JF. Statistical models and methods for lifetime data. New York: Wiley; 1982. Meeker WQ, Escobar LA. Statistical methods for reliability data. New York: Wiley; 1998. Elandt-Johnson RC, Johnson NL. Survival models and data analysis. New York: Wiley; 1982. David HA, Moeschberger ML. The theory of competing risks. London: Griffin & Company; 1978.