Contending Remaining Useful Life Algorithms

Operational readiness, reliability and safety are all enhanced through condition monitoring. That said, for many assets, there is still a need for a prognostic capability to calculate remaining useful life (RUL). RUL allows operation and maintenance personnel to better schedule assets, and logisticians to order long lead time part to help improve balance of plant/asset availability. While a number of RUL techniques have been reported, we have focused on fatigue crack growth models (as opposed to physics or deep learning of based models). This paper compares the performance of stress intensity models (linear elastic model, e.g. Paris’ Law), to Head’s theory (geomatical similarity hypothesis) and to Dislocation/Energy theories of crack growth. It will be shown that these models differ mainly in the crack growth exponent, and that this leads to large differences in the estimation of RUL during early state fault propagation, though the results of all three models converge as the RUL is shorted.


INTRODUCTION
This paper should be taken in the context of the Health and Usage Monitoring System (HUMS), which supports condition monitoring of helicopters. HUMS typically include functions for engine performance checks, automated reporting of limitations/exceedances, rotor track and balance, and aircraft regime (in order to determine when it is appropriate to acquire data).
In general, helicopters have inspections ever 50 hours of flight time, with heavier maintenance being conducted at 100 and 300 hours. Aircraft also have annual inspections. Typically, the number of hours flown per month is very dependent of the operator's mission. It is not surprising to see fleets that average 300 to 500 hours per annum. Of course, for operators that fly inspection (inventorying power poles and examining power lines for encroachment) or other missions which are seasonal (firefighting, by dropping water or delivering man/material), these aircraft can fly as much as 25 to 40 hours per week.
The importance of having a valid RUL calculation is not academic. It allows the fleet operation manager to order parts, schedule the right personnel to perform maintenance, and turn an unscheduled maintenance action into scheduled maintenance. For example, as in Figure 1, this bearing has been trending for 150 hours. The operations manager knows that in 30 hours the aircraft has a 100-hour inspection. While the aircraft is down for maintenance, the maintenance manager will schedule the bearing to be replaced. This prevents a potential unscheduled maintenance activity in the future.

Figure 1 Trending Upper Mast Ball Bearing with RUL of 83 Hours
For the light helicopter market, or cases where the aircraft has no extended overwater flights, HUMS is more of a logistic support tool to improve availability (e.g. allowing the generation of revenue flights) rather than safety of flight requirement. The aircraft, having been type certificated and properly maintained, is inherently safe. For the bearing in Figure 1, the worst possible outcome would be that the pilot sees a chips light (annunciator indicating metal debris in the gearbox) and is forced to land. Obviously, one goal of HUMS is to generate a maintenance action prior to a chips light, so that this maintenance is done opportunistically while the aircraft is down for some scheduled inspection.

RUL: INHERENT REQUIREMENTS
Inherent in calculating RUL are at least four pieces of information.
1.  (2005) gives an excellent review of many common condition indicators.
The estimate of when it is appropriate to do maintenance is a threshold setting problem. We have approached this as hypothesis testing (Bechhoefer, 2007(Bechhoefer, , 2011

Threshold Setting
All condition indicators (CIs) have a probability distribution function (PDF). Any operation on the CI to form a health index (HI) is then a function of distributions. The HI function for this paper is defined as the norm of n CIs (e.g. the normalized energy of n CIs): where Y is the whitened, normalized array of CIs, and crit, is the critical value. In hypothesis test, the critical value is calculated from the inverse cumulative distribution function, for a given probability of false alarm. For Eq. (1), the ICDF is the Nakagami where h is the number of CIs in the array, and = n, and w = h/(2-p/2)*2, see Bechhoefer et. al (2011) for the proof. From a hypothesis testing perspective, a normalized HI > 0.35 indicates a rejection of the Null Hypothesis that the component is nominal. These threshold values have been set based on the experience of monitoring numerous helicopters, wind turbines and seeded fault testing on 60+ gearboxes. The level of damage for an HI of 1.0 is seen in figure 3.
This function (1) is valid if and only if the distribution (e.g., CIs) are independent and identical (e.g., IID). As an example, for Gaussian distribution, subtracting the mean and dividing by the standard deviation will give identical Z distributions. The issue of ensuring independence is much more difficult. In general, the correlation between CIs is non-zero. This has been measured on numerous tests, see Bechhoefer et. al. 2007Bechhoefer et. al. , 2011 This correlation between CIs implies that for a given function of distributions to have a threshold that operationally meets the design probability of false alarm (PFA), the CIs must be whitened (e.g., de-correlated). In Bechhoefer et. al (2011) it was shown that a whitening solution can be found using Cholesky decomposition.
The Cholesky decomposition of Hermitian, positive definite matrix results in A = LL*, where L is a lower triangular, and L* is its conjugate transpose. By definition, the inverse covariance is positive definite Hermitian. It then follows that if: The vector CI is the correlated CIs processed as a result of data acquisition on the aircraft, which are used for the HI calculation. The transformed vector Y is 1 to n independent CIs with unit variance (one CI representing the trivial case).
The Cholesky decomposition, in effect, creates the square root of the inverse covariance. This, in turn, is analogous to dividing the CI by its standard deviation (as in the case of one CI). It can be shown that Y = L x CI T then creates the necessary independent and identical distributions required to calculate the critical values for a function of distributions.
The critical (crit, eq. 1) value is taken from the inverse cumulative distribution function for the HI. The CIs used have are assumed to have Rayleigh like PDFs (e.g., heavily tailed). For magnitude-based CIs, it can be shown that for the nominal case, the CI probability distribution function (PDF) is Rayleigh (Bechhoefer, 2005). For Gear CIs and Bearing CIs (where magnitudes which are biased by root mean square (RMS)), a transform is used to make the CI more Rayleigh like. The transform "left shifts" the CI. For example, a shift such that the .05 CDF (cumulative distribution function) is assigned to 0.0.
Consequently, the HI function is based using the Rayleigh distribution. The PDF for the Rayleigh distribution uses a single parameter, b, defining the mean µ = b*(p/2) 0.5 , and variance s 2 = (2 -p/2) * b 2 . The PDF of the Rayleigh is: x/b 2 exp(x/2b 2 ). Note that when applying these operations to the whitening process, the value for b for each CI will then be: For the HI equation in (1), the normalized energy of the CIs, it can be shown that the function defines a Nakagami PDF (Bechhoefer et. al, 2011). As note previously, the statistics for the Nakagami are h = n, and w = 1/(2-p/2)*2*n, where n is the number IID CIs used in the HI calculation.

The Concept or RUL
To be clear, for HUMS, RUL is taken as when it is appropriate to do maintenance, and not when it fails. For aviation application, maintenance is a process to restore the equipment to the original design reliability. Warn or damaged parts have reduced reliability -and as such, maintenance triggered by an HI exceeding 1, is complimentary to existing maintenance practices. For an example of a critical system, the design reliability is typically "six-nines," e.g., the probability of failure of the part under design loads is less than 10 "# per hour. For the damaged part, the reliability may be reduced to three-nines or a probability of failure of 10 "$ . The appropriateness to repair the faulty component, then, can be seen as an action to restore the designed reliability of the system as a whole. From a maintainer perspective, then: • HI ranges from 0 to a positive value, where the probability of exceeding an HI of 0.35 is the PFA.
• A warning (yellow) alert is generated when the HI is greater than or equal to 0.75. Maintenance should be planned by estimating the RUL until the HI is 1.0.
• An alarm (red) alert is generated when the HI is greater than or equal to 1.0. Continued operations could cause collateral damage. This threshold setting model ensures that the probably of false alarm is exceedingly small when the HI reaches 1. In practice from numerous installations and seeded fault test, a bearing at HI 1 has easily seen physical damage (figure 3).
Note that this nomenclature does not define a probability of failure for the component, or that the component fails when the HI is 1.0. Instead, it suggests a change in operator behavior to a proactive maintenance policy: perform maintenance before the generations of collateral or cascading faults. By performing maintenance on a bearing prior the bearing shedding extensive material, costly gearbox replacement can be avoided, and the reliability of the gearbox is restored to its design requirements.
Thus, we define the RUL as the time from the current HI, until the HI is greater than or equal to 1.

Estimating Future Loads
If the cyclic loading/stress level is sufficiently high, microcracking will propagate across the surface of the component, and eventually penetrate into the body of the material. For any material, to predict the RUL, an estimate of the future load is needed. If the future load is below the plain fatigue limit of the material, the component will run indefinitely.
For HUMS, there is limited knowledge of load except for gross measure such as torque (Figure 2), or airspeed. For many accessory components (oil pump, hydraulics), the loads are relatively unchanged by load or mission. What can be said is that some missions require more power/load than other missions, based on factors such as altitude, temperature and regime (hover requires more power than level flight at 40 knots, etc.).

Figure 2 PDF of Torque on a Light Helicopter
Because of the limited knowledge available for strain, it is proposed that three cases would be used to bound the RUL corresponding to: low, medium and high loads. The low estimate is 10% lower than the medium load/torque (e.g. average mission), whereas high load is 10% higher.
It was found that the mean load was easily calculated after 20 flight hours.

MODELING COMPONENT DEGRADATION
Generally speaking, fatigue crack growth can be characterized as (Beer, 1992): Mode1 (opening mode), the crack surface moves directly apart, Mode 2 (edge sliding mode), where the crack surfaces move normal to the crack front, remaining in plane, and, Mode 3 (shear mode), the crack surface moves parallel to the crack front and remains in the crack plane.
These three modes can describe, through superposition, most general cases of surface displacement. It is assumed that cracks are stress-free boundaries adjacent to the crack tip. As such, loading forces affect only the intensity of the stress field at the crack tip. These fields correspond to the three modes of crack surface displacement and are characterized by the stress intensity factor K (which is a function of the component dimension and loading condition). K is proportional to the gross stress and the square root of the crack length. The opening mode intensity factor is given by: Where s is the gross stress, a is the crack length and a is the shape factor. When K is known, stress and displacements near the crack tip can be calculated. For example, for a displacement in the x direction, the strain is (see Frost, 1999): This representation of crack tip stress field by a stress intensity factor is a basic concept in fracture mechanics. In the analysis of fatigue crack growth, the fatigue cycle is usually described by DK, which is the difference between the maximum and minimum stress field. It has been seen, experimentally, that DK, and not the maximum K, causes crack growth, and that if DK is constant, the crack growth rate is constant (Frost, 1999).

A Linear elastic fracture mechanics model
For many materials, such as steel used in gears and bearings, which are subject to tensile loading cycle, the fatigue crack growth can be expressed as: where • da/dN is the rate of change in the half crack length per cycle

• D is a material constant
• m is the crack growth exponent, typically 3 to 5.
Substituting in DK : Inverting and integrating to get N, the number of cycles gives: By taking a as ao to get the crack growth rate, the constants cancel out leaving: If m is allowed to be 4, this give: For constant rate machines, such as a helicopter gearbox, N is proportional to time.
If one makes the assumption that the component health (the HI) is proportional to damage, the Eq. (9) defines the RUL estimate.

Head's Theory as a fracture mechanic model
Head's theory makes a number of simplifying assumptions about the stress field. Essentially, material near the crack Is treated as an array of independent elastic tensile bars of modulus E, each carrying the remotely applied stress s, which transmitted the load to the bars both directly and through shear. The model gave for an applied stress s, This is interesting because in effect, other than a constant, Eq. (10) can be reduced to Eq. (6), where the exponent m is 6. Inverting and taking the integral gives the following alternate RUL algorithm:

Dissociation theory as a potential fracture mechanics model.
In the case where crack loading is in the anti-plane strain (e.g. Mode 3), the plastic zone at the crack tip can be represented as a continuously distributed array of small dislocations at on the crack plane. It is assumed that crack growth occurs when the accumulated plastic strain distribution at the crack tip exceeds some critical value and continues as this value is exceeded at the crack tip. The rate at which the crack grows per stress cycle in terms of displacement leads to: This again is similar to Eq. (6), with an exponent of 2. Inverting, integrating and changing terms gives:

COMPARISON OF THE THREE FRACTURE MECHANICS MODELS
In order to compare the performance of these three models, a data set is needed. We used data acquired from a 2.2 Megawatt (MW) high speed bearing, with a known fault (HI value 1.1, Figure 3). This data set was collected over 55 days, one acquisition every 10 minutes (144 acquisitions per day).
Note that the fault starts to propagate at approximately time -700, which corresponded to high loads from a winter storm ( Figure 4). All three of the fracture mechanics models require a crack length (a), and the rate of change of the crack length (da/dN).
A simple state reconstruction model (α-β tracker as in "alphabeta") was used to filter the component health (HI) and estimate the rate of change of the component health (e.g. dHI/dt). Figure 5 shows the estimated rate of change of health using the α-β tracker.
By making the simplifying assumption of stationarity, we will substitute the α-β tracker for Kalman filters or some other Ricote equation.

Figure 4 Measured and Filtered HI for a High-Speed Bearing
As such, taking the limit as time moves toward infinity of the Kalman filter, the coefficients for the α-β tracker (used for HI, and dHI/dt) can be calculated as a constant (Bar-Shalom, 1992): Where the process variance is sw 2 , and plant noise variance is sv 2 . The filter gains are: The state equation are then: For each HI update: fHI = fHI + dHI * dt; rk = HI -fHI; fHI = fHI + alpha*rk; dHI = dHI + (beta*rk)/dt; Observer in figure 4 and Figure 5. that before -700 hours, the HI, and DHIT/dt is relatively flat, and in effect the RUL would be large or infinite. One would expect that as the fault starts to propagate, the RUL reduces quickly, as can be seen in figure 6.

Comparing the Models
It is assumed in this construct that there is no physical measure of component damage, e.g. crack length. It is proposed that the current crack length is proportional to the current HI, with that measurement being corrupted by additive Gaussian noise. Further, the crack rate of growth, is dHI/dt. Using this definition, the RUL values for the three models, derived from Eq. (9), (11), and (13), are seen in Figure 6.

Figure 6 Comparison of RUL models
The ideal model is one in which the RUL is the number of hours prior to perform maintenance. The slope of the RUL is -1 hr/dt, e.g. one hour of life is consumed for each are the machine is run. This an idealized construct is based on an average future load. Note that from the start of the data acquisition until time -700, the rate change of health if effectively zero. As such a minimum degradation rate is used (1 over the mean time between failure). This is done numerically to prevent a divide by zero error. Figure 7 Shows a model performance from the start of the fault propagation until the HI value of 1 is reached.
As observed, all three models converge and would provide useful information to a maintainer or aircraft scheduler. Conceptually, the model with the best dynamics (best fit) should have a derivative that approaches the idealized RUL model of -1 hr/dt. Using an α-β tracker, the derivate of each model was calculated to evaluate the model dynamics performance (Figure 8). Both linear and dislocation models seem to outperform Head Theory model. In term of quantifying performance, we us the simple mean, median and RMS of the derivate from say, time -350 to 0 (Table 1). It is clear that modeling errors exists. Sources of error maybe that that the crack growth exponent may be a non-integer value between the linear elastic model and dislocation theory model. Alternatively, the CIs may not be linearly related to crack growth. Note that for an ideal RUL model, the mean and median for the derivative of RUL would be -1. As the linear elastic model has better performance than either the Heads or Dislocation model, it can be stated that fatigue crack growth can be characterized as Mode1 (opening mode), where the crack surface moves directly apart.

Predicting the Future Component Health
One can compare the predicted component health by solving for af for some given time in the future. For example, from equation (9), the predicted HI for the linear elastic model would be: Where da/dn is the current calculate HI derivative and ao is the current calculate filtered HI (Figure 9). Figure 9 shows an example at time -400 hours (that is, given that in 420 hours, the filtered HI is at 1, as we have this knowledge post processing), we can see that the actual trajectory of the fault propagation lies somewhere between the linear elastic (RUL of 429, or 7% error), and the dislocation model (RUL of 360 hours or -10% error), whereas Heads theory has an RUL of 463 hours (16% error). Figure 9 Predicted HI trend for the three models

Estimating the crack growth exponent
If one assumes from the data that at time -700 until zero, there is a fault propagation and that the RUL decreases linearly, then by rearranging terms in eq 8, one gets: In an ideal linear RUL model, one could solve for m with (18) using the calculated dHI/dt and HI values with a0 equal to the current HI, with af set to 1. It would be expected that the exponent is between 2 and 6, given the results from Table 1, however, the results of this were mixed ( Figure 10).

Figure 10 Calculated exponent from data
Conceptually, the rate of change of health and any time index t (e.g. when we acquire data), is going to be dependent on the accumulated load (strain) from t-1 until t. This typically is not well captured when measuring data periodically and is a source of error. Other possible sources of error, or why the model exponent was not exactly 2, may be the result of lag resulting in the calculation of the filtered health, or errors in the estimation of the rate of change in health. These error sources are a continuing research interest.

DISCUSSIONS
For many owners and operators, a helicopter is a utility vehicle used to accomplish a mission. These missions' requirements are not easily met by other types of vehicles. The operator only generates revenue when the aircraft is flying. The only time a helicopter can't fly, other than weather conditions, is when it's in maintenance. Most of the maintenance activities can be scheduled. Manufacturers provide planned maintenance to restore the aircraft to its design reliability. Maintenance is done on the aircraft to replace parts (for those affected by a retirement life). Inspections, with their appropriate intervals, allows operators to find and replace wear items on the helicopter.
These are easily manageable and are scheduled in advance by maintainers. Unscheduled events and incidents (chips light) remove the aircraft from service and negate the opportunity to generate revenue. For example, a bearing that starts to wear is usually found during a scheduled inspection or overhaul of the components. If the bearing is damaged between an inspections or overhauls, the helicopter is designed to identify the fault (e.g., a "Chip" light would illuminate in the cockpit identifying an issue has occurred). These types of defects result in removing the aircraft from service to perform maintenance.
Implementation of a prognostics system on the helicopter allows the maintainer/operator to identify the component, converting an unscheduled maintenance event into scheduled maintenance. A prognostic system, giving a maintainer an estimate of the Remaining Useful Life (RUL) allows the following: · Operations to schedule the aircraft for maintenance opportunistically with a current event.
· Maintainer to order parts to marshal any specialized tool/capabilities needed for the repair.
For helicopter manufacturers supporting customers around the world, the goal is to keep customers helicopter in the air ensuring it is safe to fly. One of the hardest things to do is to say to an operator, "we do not recommend that you continue to fly, this part should be replaced before next flight." This is especially the case when the aircraft is needed for a search and rescue or some other critical mission: a recommendation to not fly means that the aircraft is not mission capable.
The Prognostics health monitoring system should provide specificity of the damaged component and an estimate of RUL. Giving owners, operators, and maintainers these types of tools, allows for improved operational readiness, more opportunities for revenue, and an enhanced sense of safety.
As noted, a RUL allows the maintainer to determine which tools, capabilities, and parts are needed for a repair. Unscheduled events can cause unanticipated aircraft grounding if spare parts are not readily available. Many operators fly on a contract with a penalty clause. Unless they have spares in their inventory, they will need to find the component from the OEM or another source, place the order and have the part shipped to you. The cost of replacement components is very dependent on availability and urgency.
Having the availability to see the component wear before it starts failing is a huge advantage, improving logistics and reducing spares. Knowing the RUL of the component gives you time to order your part in advance, which reduces the cost and time to schedule the replacement.

CONCLUSION
The use of a fracture analysis model gives insight into the development of a remaining useful life model. In this paper, three high cycle fatigue models where tested and compared with respect to their ability to estimate the remaining useful life of a bearing. These models were used to measured vibration data which was collected over 55 days on a wind turbine. The vibration data was used to extract features that are correlated with bearing damage. These condition indicators were fused into a health indicator. The health indicator was constructed such that when the HI is 1, it is appropriate to do maintenance. The RUL was then the time from the current HI until the HI is 1.
All three models (linear elastic fracture mechanics, Head's theorem and a dislocation theorem model) gave satisfactory result, the greatest difference being the RUL estimate prior to propagation of the fault.
It was hypothesis that the "best" model would have a derivative of the RUL of -1. Even though the degradation process is nonlinear, the RUL, based on some mean load, is linear: In a high confidence model, for each hour consumed in operation, the RUL would decrease by 1. Even in the highly variable loads seen in a wind turbine, it was shown that is a valid assumption.
It was found that the linear elastic model has a mean, and median derivative close to one, confirming that it is the best model. That is, the fatigue crack growth can be characterized as Mode1 (opening mode), the crack surface moves directly apart The Heads theorem RUL derivative was greater than one (indicating it was overestimating the RUL). The dislocation theorem model RUL derivative was less than 1, indicating that it was underestimating the RUL.
As the linear elastic model is shown to be the most accurate, then the crack growth exponent value should be approximately 2. Calculating this from the data, the exponent ranged between 3 and 1.
All three models could detect the fault propagation and gave an acceptable RUL (e.g. meaning the data could be used for maintenance planning) 400 hours prior to when it was appropriate to do maintenance.
condition-based maintenance, rotor track and balance, vibration analysis of rotating machinery and fault detection in electronic systems. Dr. Bechhoefer is a board member of the Prognostics Health Management Society and a member of the IEEE Reliability Society.

Marc Dubé received his College Degree in Aircraft
Maintenance from the École Nationale D'Aérotechnique (ÉNA). He is a TCCA licensed Aircraft Maintenance Engineer who has several years of experience in helicopter maintenance. From field maintenance to maintenance management roles with different Canadian operators/maintainers, he has been part of the Bell product Support Engineering (PSE) team for a number of years. Specializing in Rotors, Track and Balance, Drivetrain and HUMS, Marc supports Bell customers around the world.