Determination of the Optimal Replacement Age for the Preventive Maintenance of Bearing Assemblies Involving Weibull Failure Probability Distribution Functions

This paper describes the determination of optimal replacement age for bearing systems with failure times described by Weibull distributions. Optimal age replacement policies are determined for bearing assemblies produced two different types of steel. The parameters of the Weibull probability distribution functions for the bearings were determined by non-linear least square analysis from published data on rolling contact fatigue lives. The resulting distribution functions are used as inputs into the standard expression for the maintenance cost of an age replacement policy and manipulated symbolically using the computer program Maple. These yields closed form expressions for the policy costs that invariably exhibit the well-known vase shape characteristic of these types of problems. The resulting expressions can then be easily used to determine the optimal replacement age of the bearing components. The problem of determining the optimal age replacement policy of bearing assembles consisting of independent components arranged in series is also examined. The effect on the optimal age replacement time of using the same type steel to manufacture all the components is compared with that of building the bearing assembly using components made with different steels. As expected, the results clearly show the increased superiority of the higher-quality steel components in the form of much longer optimal replacement ages. However, replacement policy efficiencies depend on replacement time in a complex fashion. Moreover, the results also suggest that building bearing assemblies combining high quality and low quality steel components may be wasteful. Overall, computer experimentation and examination of the behavior of cost functions using symbolic manipulation with Maple can produce useful guidelines for the design of optimal age replacement policies.

Abstract This paper describes the determination of optimal replacement age for bearing systems with failure times described by Weibull distributions. Optimal age replacement policies are determined for bearing assemblies produced two different types of steel. The parameters of the Weibull probability distribution functions for the bearings were determined by non-linear least square analysis from published data on rolling contact fatigue lives. The resulting distribution functions are used as inputs into the standard expression for the maintenance cost of an age replacement policy and manipulated symbolically using the computer program Maple. These yields closed form expressions for the policy costs that invariably exhibit the well-known vase shape characteristic of these types of problems. The resulting expressions can then be easily used to determine the optimal replacement age of the bearing components. The problem of determining the optimal age replacement policy of bearing assembles consisting of independent components arranged in series is also examined. The effect on the optimal age replacement time of using the same type steel to manufacture all the components is compared with that of building the bearing assembly using components made with different steels. As expected, the results clearly show the increased superiority of the higher-quality steel components in the form of much longer optimal replacement ages. However, replacement

Introduction
Striving for savings and efficiency, industries often sacrifice maintenance. However, the neglect of proper maintenance of complex technological systems can involve significant risk. As an example, it has been reported that over 20% of medical device malfunctions are due to incorrect maintenance (Servomaa, 1989). Also, almost 12% of all aircraft accident reports cite maintenance as a factor and between 1994 and 2004; maintenance problems have contributed to up to 42% of fatal airline accidents in 288 Determination of the Optimal Replacement Age for the Preventive Maintenance of Bearing Assemblies Involving Weibull Failure Probability Distribution Functions the United States alone (O'Brien 2012). Three common types of maintenance are preventive, corrective and failure-finding. In preventive maintenance one attempts to reduce the failure probability of the system by acting even when the item is functioning properly. Preventive maintenance can be age-based, clock-based, condition-based or opportunity based. In age based maintenance, replacement takes place once the system reaches a predetermined operating age. Clock based maintenance is applied after a predetermined length of time. Condition based maintenance is applied when continuously measured diagnostic signals suggest that the system needs replacement, regardless of its age or clock time. Corrective maintenance is really repair and takes place once the item has failed and it may be impossible in mission critical systems. Failure-finding maintenance is the result of failure detection during continuous monitoring. Again, this cannot be afforded in mission critical systems.
The optimum age replacement policy is the one that minimizes the average expected cost per unit time over an infinite time span. This paper describes a simple method for the determination of the optimal replacement age for preventive maintenance of a system whose failure probability is described by a Weibull distribution. Two and three parameter Weibull distributions were investigated. The symbolic manipulation program Maple (Maplesoft) was used to obtain solutions of the problem. For certain values of the parameters in the distribution, convenient closed form expressions are obtained. Numerical methods must be applied to obtain solutions for general values of the parameters.
Preventive age maintenance is often used as maintenance policy for the proper operation of bearing assembly systems. Under this policy, if an assembly survives to a specified age t 0 , it is replaced and if the assembly fails prior to attaining the specified age t 0 , it is replaced and the new unit is scheduled for replacement if and when it reaches age t 0 . Clearly, replacing bearing assemblies too frequently is costly and wasteful since most bearings would still be in good condition while waiting too long before replacing, risks the occurrence of catastrophic failure, which may be quite costly. Hence, the problem of determining the optimal replacement age has long attracted a great deal of interest.
The study of optimum preventive maintenance policies has long been the subject of theoretical and practical interest (Lotka 1939, Campbell 1941, Barlow and Proschan 1964, 1965, Fox 1966, Glasser 1967, Subramanian and Wolff 1976, Berg 1976, Nakagawa 1977, 2005, Bergman 1982, 1985, Valdez-Flores, C. and Feldman, R.M. 1989, Sandve and Aven 1999, Sheu, S-H, and Griffith W.S. 2001. Coolen-Schrijner, P. and Coolen, F.P.A. 2006, Mahdavi, M. and Mahdavi, M. 2009, Sarkar et al. 2011). This work further explores age replacement maintenance policies by focusing on the use of symbolic manipulation to investigate the specific problem steel ball bearings and obtain closed from expressions of the policy cost as a function of the replacement age, the values of the parameters in the failure probability distribution function and the potential liability costs associated with unexpected failure. To the author's knowledge, this is the first instance of the use of symbolic manipulation techniques in the investigation of age replacement policies.

Description of the System Investigated
Rolling bearings are machine elements used to constrain motion in a desired direction while carrying load under low friction conditions between contacting parts in relative movement (Hamrock and Anderson 1983). Bearings have been used in practical applications for thousands of years and played a key role during the Industrial Revolution (ABMA undated, Dowson 1998). Bearing design has evolved over centuries and bearings are nowadays sophisticated indispensable elements allowing the efficient operation of modern industrial machinery.
Many different types of bearing are encountered in applications. Journal, ball and roller bearings are widely used. All bearings are really assemblies consisting of several components. For example, figure 1 is a photograph of a partially disassembled roller bearing showing the steel inner and outer rings (races), steel rolling elements and brass cage.
During their normal operation, the various bearing components, especially the rolling elements and race surfaces, experience the application of high stresses repeatedly over small volumes of material. The cyclic loading over long periods of time eventually may result in the dislodging of pieces of material by a process known as surface fatigue failure. Failure of any of the components hampers bearing performance and often results in failure of the assembly.   Figure 2 shows an example of bearing assembly failure due to the formation of a scar on the surface of a journal bearing sleeve, most likely produced by a fatigue crack originating near the surface of the sleeve that grew large enough that a flake-like piece of material was eventually removed (Desjardins 2014). The scars resulting from flaking or spalling of material from the surfaces of bearings produced undesirable features in the pressure distribution and operating parameters of the bearing assembly and lead, in the best cases, to unwanted stoppages and the need for corrective maintenance and in the worst case, to potentially catastrophic outcomes.

Statistics of Ball Bearing Failure
It is generally accepted that failure times of bearing assemblies and bearing assembly components are often well represented by Weibull distributions (Lieblein andZelen 1956, Zaretsky andPoplawski 2000). Two and three parameter Weibull distributions were introduced (Weibull 1939(Weibull , 1951 who demonstrated their flexibility and wide applicability. The Weibull distribution was soon used as basis to obtain now classical expressions for the prediction of bearing assemblies lives (Lundberg and Palmgren 1947;Gupta et al. 2013) The two parameter Weibull failure probability distribution function is defined by the expression Where α and β are the location and shape parameters of the distribution. This function is flexible and of general applicability in reliability analysis. The parameters α and β are always positive numbers and F(0)=0. Moreover, in the limit of β=1 it becomes the exponential distribution and as β increases, it resembles first a normal distribution, then a Log-normal distribution with long tail.
The survival distribution or reliability function is the complement of the failure distribution and is defined as S = 1 -F Survival of the bearing assembly requires survival of all the components and since the failure of one component occurs independently of the failure of another, the probability of survival of the assembly can be estimated as the product of the survival probabilities of the individual components (Mayer 1970), e.g. for a ball bearing assembly consisting of three components, inner race, outer race and balls, the survival function of the assembly is S = S ir S or S b Here, S ir , S or and S b are the survival functions of inner race, outer race and bearing balls, respectively. Even if the independence assumption is invalid, the failure times of bearing assemblies are often also reasonably well described by a Weibull distribution. It is herein then assumed that the failure probability distribution function F of a bearing assembly or any of its components is reasonably well represented by a two parameter Weibull distribution.
A nice feature of the Weibull distribution is that simple closed-form expressions are readily obtained for important reliability metrics such as the probability density function f, given by Where Γ(1/β) is the Gamma function of argument 1/β (Evans et al. 2000).
In short, much is known about the reliability of a specific type of bearing assembly once one has determined the values of the parameters α and β for that particular bearing assembly.

Formulation of the Optimum Preventive Age Maintenance Problem
Following is the formulation of the optimum age replacement problem for a bearing assembly.
If the assembly or any of its components is replaced before it fails, the replacement cost incurred is c. On the other hand, if one waits until the item fails before it is replaced, in addition to the replacement cost c one must also pay extra expenses (e.g. liability costs) k that, depending on the criticality of the component, can be several times larger than the replacement cost (i.e. k = r c where r >> 1 is an empirical proportionality coefficient).
The goal is to determine the optimal replacement age t 0 * of the system under an age-based preventive maintenance policy. Let T be the time to failure of the bearing assembly and t 0 be the replacement age of the assembly (i.e. the time at which the assembly is to be replaced). As mentioned above, the survival probability of the system is assumed to be well represented by the product of the survival probabilities of the constituent components.
If T is the time to failure of the bearing assembly in millions of cycles, the mean time between replacements of assemblies with replacement age t 0 is given by (Rausand, M. and Hoyland, A. 2004, Sec. 9.6) where f(t) is the failure probability density function, Pr(T >t 0 ) is the probability that failure would not take place before replacement and F(t) is the failure probability distribution function. The total cost per replacement period is equal to the replacement cost plus the extra cost in case failure takes place, i.e. c + k F(t 0 ) and the total mean cost per unit time ($/million cycles) is equal to the total cost per replacement period divided by the mean time between replacements at age t 0 , This is exactly the expected cost per unit time in the long run obtained from a fundamental theorem of renewal theory. If one has determined the specific mathematical form and parameters of the failure probability distribution function F(t), all one needs to estimate the value of the expected cost per unit time for any replacement age t 0 are the values of the replacement cost c and the multiplier coefficient r. If only the value of the replacement cost is available, (i.e. the cost of replacement parts is readily available while the costs incurred if failure occurs can be more difficult to estimate) the cost C can be represented as a function of two independent variables, t 0 and k (or r).
Note that as t 0 goes to infinity (i.e. replacement only after failure -corrective replacement), the total cost per replacement period ) dt is the mean time between failures. Clearly, if k is large then C ∞ will be large too. Note also that, as expected, when t 0 goes to zero (frequent, premature replacement), the cost of replacement becomes directly proportional to the number of times the assembly is replaced and it can become large too. Therefore, one may expect that there will be a certain item age t 0 * that would yield the lowest replacement cost per replacement period (Rausand and Hoyland 2004, Nakagawa 2005, Osaki 2002. A useful metric to judge the appropriateness of an age replacement policy is the cost efficiency of the policy E p ; this is defined as E p = C(t 0 )/ C ∞ Low values of E p indicate high cost efficiency of the policy.

Computational Experiments and Results
The software Maple (Maplesoft.com) was used to perform computer experiments through symbolic manipulation and obtain exact solutions to a couple of selected problems. Empirical data about rolling contact fatigue lives (in millions of cycles) of rolling bearings made from AISI 52100 steel by air melting (AM) and by the vacuum arc refining-electroflux remelting (VAR/EFR) process was used to define useful ranges of interest for computer experimentation (Oakes 1981). The range of lives thus determined was between 0 and 100 million cycles Non-linear least squares analysis of the published experimental data was used to obtain values of the Weibull distribution parameters a and b for both the AM and the VAR rolling bearings. The results are shown on the following table. Table 1. Computed values of location and shape parameters of the Weibull distribution for AM and VAR rolling bearings from the data in (Oakes 1981 3a,b,c,d show, respectively plots of the functions f, F, S and h for the steel bearings manufactured from both types of steels and clearly demonstrate the significant increase in durability of the bearings produced by the VAR process.      The production of VAR steel is a more complex steelmaking process than the production of AM steel. Correspondingly, the price of VAR steel can be up to ten times higher than that of AM steel. We assume that the bearing manufacturing costs for both types of steel are similar but take into account the price differential by assuming that the price of VAR steel is ten times larger than that of AM steel (i.e. c var = 10 c am ).
The Maple software can be used to determine the specific functional form of the cost functions for both types of bearings. For instance, if the value of c am = c 1 =10 for the 292 Determination of the Optimal Replacement Age for the Preventive Maintenance of Bearing Assemblies Involving Weibull Failure Probability Distribution Functions AM steel bearings and c=100 for the VAR bearings are used, the resulting cost functions C(t 0 ,r) are; for the AM bearings (C AM = C 1 ) And for the VAR bearings (C VAR = C 2 ) Figure 4a shows the shape of the computed C(t0,k) surface for the AM bearings and figure 4b shows the result for the VAR bearings.  Table 1 and c=10 (AM bearings).  Table 1 and c=100 (VAR bearings).
While the two surfaces are quite different from each other they have several common features. In both cases, the cost goes up quickly the more premature the replacement. As the replacement age increases for a fixed value of r, the cost decreases and the functions show a valley, although the valley is bigger and shallower in the case of the AM bearings.
For fixed, large values of t 0 , the cost increases rapidly with the value of r as the liability charges become greater. Figure 5a shows the computed values of the cost as a function of t 0 for values of k = 10 2 , 10 3 and 10 4 for the AM bearings and Figure 5b shows the results for the VAR bearings.  Table 1 and c=10 (AM bearings).

Figure 5b.
Cost as a function of t 0 for values of k = 10 2 , 10 3 and 10 4 , using the values of α and β in Table 1 and c=100 VAR bearings).
As expected, in both instances, when the value of the liability cost is small, optimal replacement times are larger in both cases but more so in the case of the VAR bearings i.e. replacement can be delayed without much cost penalty when using the higher quality steel. On the other hand, when liability costs are high, the cost of postponing replacement beyond the optimal replacement time increases rapidly in both cases.
The cost function for the AM bearings has relatively small values for t 0 < 6x10 6 but it increases very rapidly as t 0 increases especially for large values of k. The cost per unit time is certainly larger for the VAR bearings at small values of t 0 but it increases slowly as t 0 increases. In both cases, from a cost point of view, there is strong motivation for replacement at age around the optimal replacement age in order to both, reduce the risk exposure due to failure and to avoid wasteful replacement of perfectly functioning components.
The optimal replacement time t 0 * for any given value of k is obtained is the value of t 0 at which the extremum condition is satisfied, i.e. dC/dt 0 (t 0 *)= 0 Figures 6a and 6b shows the computed values of the derivative of the cost with respect to the replacement time as a function of replacement time for the two types of bearings for values of k = 10 3 ,10 4 ,10 5 ,10 6 . IN the case of the AM bearings, he optimal replacement time decreases from ~ 5x10 6 to about 2.7x10 6 as k increases in the range considered, whereas in the case of the VAR bearings the optimal replacement time is much more sensitive to the increase in k and decreases from 16.7x10 6 to about 0.3x10 6 .  Table 1 and c=10 (AM bearings).

Figure 6b.
Derivative of the Cost with respect to replacement time as a function of t 0 for values of k = 10 3 ,10 4 ,10 5 , 10 6 , using the values of α and β in Table 1 and c=100 (VAR bearings).
For low and moderate values of liability costs the VAR bearings are much more forgiving and allow much longer lives without replacement. For instance, when r=10 the optimal replacement time rises to about 15x10 6 cycles for the VAR steel while it only goes up only to about 6.5x10 6 for the AM bearing. However, in case the liability costs are high, optimal replacement times are short regardless of the quality of the materials involved (4-5 x10 6 cycles). While the same value of replacement cost c=10 has been assumed for components made of the two steels, this is clearly not the case in practice since VAR steel is generally more 294 Determination of the Optimal Replacement Age for the Preventive Maintenance of Bearing Assemblies Involving Weibull Failure Probability Distribution Functions costly. However, the computed value of t 0 * for the VAR steel bearings changes only slightly even if the value of c for the VAR steel triples. If the more costly VAR steel bearing is used, the higher replacement cost involved may make its use disadvantageous against the conventional AM bearings when the liability cost is high. In such a case, the liability risk overshadows the benefit of using the higher quality material and early replacement using the least expensive component seems to be the safest and cheapest strategy.
The cost efficiencies of age maintenance policies for bearings manufactured from AM and VAR steels for k=10 6 are shown in figure 7. The figure shows that both policy efficiencies decrease as t 0 increases. However, for values of t 0 less than about 8x10 6 , the age replacement policy for AM steel bearings is more efficient than the one for VAR steel bearings. However, as t 0 increases beyond 8x10 6 cycles, the efficiency of the policy when using AM steel bearings decreases rapidly and becomes less efficient than the age replacement policy when using the VAR steel bearings. In other words, if replacement ages are chosen below 8x10 6 cycles it is most cost efficient to use AM steel bearings, otherwise, use of the VAR steel bearings is more cost efficient.. Note that the above expressions for the cost functions give cost values explicitly in terms of the values of the independent variables of the problem (t 0 and k). However, the expressions are cumbersome and invoke special functions which may make their use without the benefit of symbolic manipulation with Maple tedious and somewhat impractical for in everyday reliability applications. In the sequel, an approach to produce solutions of simple form is presented that can remedy the above drawback to some extent.
As another illustration of the use of symbolic manipulation in preventive maintenance analysis, consider the case of a simple journal bearing assembly but now assuming that assembly failure can be due to failure of any of the two main components (journal/shaft and bearing/sleeve). Under the assumption that failure of any component represents failure of the whole assembly, the failure time probability distribution function of the bearing assembly is then given by Where α j , β j and α b , β b are the Weibull distribution parameters for the journal and the bearing, respectively. The distribution function F for the journal-bearing assembly in which the two components are made of VAR steel follows a Weibull distribution with parameters α = α var /(2 (1/βvar) ) = 32.2640and β = β var = 1.7896, i.e. F = 1 -exp(-(t/32.2640) 1.7896 ) The symbolic manipulation software takes a long time attempting to produce a closed form expression for the cost. However, it is possible to quickly obtain simple expressions that are good estimates of the cost by using rounded, integer values for the parameters. Specifically, consider the expression F = 1 -exp(-(t/32) 2 ) Figure 8 shows a comparison of the two functions above; although clearly not exactly the same, the two functions are very similar. Introducing the second expression into the cost equation, with c=100, yields in Maple the following simple closed form expression Figure 9 shows the computed cost as a function of t 0 and k for this case. Differentiation with respect to t0 yields the following closed form expression Figure 10 shows the shape of the dC/dt 0 surface as a function of t 0 and k. Figure 10. Computed values of the function dC/dt 0 versus t 0 and k assuming F = 1 -exp(-(t/32) 2 ) for c=100.
From the above, the optimal age replacement time t 0 * for k=10 4 is determined as 3.2x10 6 , which is somewhat smaller than the value of 4.2x10 6 obtained for the bearing considered as a single component. The mean time between failures is calculated as t MTBF = ∫ 0 ∞ 1 -F(t) ) dt = 28.35 million cycles The calculated value of t MTBF for the bearing system regarded as a single unit is 42.24x10 6 cycles (see figure 3e above). The smaller calculated values of t 0 * and t MTBF for the composite system are directly related to the reduction in the location parameter value of the Weibull distribution of such system resulting from the reduced reliability produced by the assumption of the two components being connected in series. Values of F versus t were computed for three cases. In the first case, both components are manufactured of VAR steel, in the second case, both components are produced using AM steel and in the third case, one component is made of VAR steel and the other of AM steel. For a given lifetime, values of the failure distribution function of the assembly are higher than those of the components. Results show that significant improvement in bearing journal-bearing assembly life can be obtained when switching from AM steel to VAR steel. However, both components of the bearing assembly must be manufactured using VAR steel because combining steels produces only a marginal increase in optimal replacement time over the use of AM steel for both components and is therefore wasteful.
A useful finding about the symbolic manipulation software is that it can produce relatively simple closed form expressions for the cost function if integer values of the location and shape parameters of the Weibull distribution. Therefore rounding the values obtained from nonlinear least squares processing of empirical data to the closest integer can be recommended as a means to produce reasonable approximations to the cost functions of actual systems with simple mathematical form.
The results obtained for all the three two component bearing assembly cases above show again clearly the rapid rise in cost obtained with too frequent replacement of young components in good condition. The results also show that the costs rise with replacement age past the minimum but then settle to the asymptotic maximum value of C ∞ in each case. In all situations investigated in this study, computer experimentation shows that the cost always appears to be smooth function of the independent variables and it seems there is always a certain component age t 0 *(r) that yields the lowest cost. Also, the optimal replacement age shifts to lower values as the liability costs increase. Even if no closed-form expressions for the optimal replacement time can be found by Maple, the software can still be used to estimate it by examining the behavior of the derivative dC/dt 0 .

Conclusions
This paper shows that solutions to optimal problems of age-based preventive maintenance where the failure probability distribution function is of Weibull type can be 296 Determination of the Optimal Replacement Age for the Preventive Maintenance of Bearing Assemblies Involving Weibull Failure Probability Distribution Functions easily obtained by computer aided symbolic manipulation tools (e.g. Maple). The approach is demonstrated by application to the life of bearing assemblies. Such assemblies can be regarded as composed of any number of components arranged in series, each with its own failure probability distribution function. In all cases examined, the cost rate (replacement cost per unit time) is obtained as a monotonic function of the replacement time t 0 and the liability cost of premature failure k. For any given value of k, the cost function exhibits a minimum value for a specific value of t 0 = t 0 *. Symbolic manipulation is capable of producing simple closed form expressions for the cost per unit time as a function of replacement time and the liability cost. Examination of the resulting cost functions can then be used to produce accurate estimates of optimal replacement ages as a function of the cost of components and the liability cost of premature failure. The methodology discussed should be applicable to the problem of optimal condition-based maintenance.