Review and application of Rainflow residue processing techniques for accurate fatigue damage estimation

Most fatigue loaded structural components are subjected to variable amplitude loads which must be processed into a form that is compatible with design life calculations. Rainflow counting allows individual stress cycles to be identified where they form a closed stress‐strain hysteresis loop within a random signal, but inevitably leaves a residue of open data points which must be post-processed. Comparison is made between conventional methods of processing the residue data points, which may be non-conservative, and a more versatile method, presented by Amzallag et al. (1994), which allows transition cycles to be processed accurately. This paper presents an analytical proof of the method presented by Amzallag et al. The impact of residue processing on fatigue calculations is demonstrated through the application and comparison of the different techniques in two case studies using long term, high resolution data sets. The most significance is found when the load process results in a slowly varying mean stress which is not fully accounted for by traditional Rainflow counting methods. 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license


Introduction
The calculation of conservative load cycle spectra is a fundamental aspect of fatigue design, requiring an estimate to be made of expected operational loading conditions. Complex lifecycle loading may be simplified by dividing the process into discrete load cases, such as take-off and steady flight conditions for the analysis of aircraft components. Fatigue life can then be quantified in terms of time to crack initiation through the concept of linear damage accumulation, or by the application of crack growth models. Both approaches utilise information about the range, mean and number of stress cycles that will occur [1].
The identification of individual fatigue loading cycles within a random stress amplitude time series is achieved through the use of a suitable cycle counting algorithm. Typical methods include level-crossing counting, range-pair counting, reservoir counting, and Rainflow counting. Variations of these algorithms are included in the ASTM cycle counting standard [2].

Background to the Rainflow counting algorithm
Rainflow (RF) counting has become the most widely accepted method for the processing of random signals for fatigue analysis, and testing has demonstrated good agreement with measured fatigue lives when compared to other counting algorithms [3]. The concept was first developed by Matsuishi and Endo [4], where the identification of cycles was likened to the path taken by rain running down a pagoda roof. In the paper, the authors defined a full RF cycle as a stress range formed by two points which are bounded within adjacent points of higher and lower magnitude; as the stress path returns past the first turning point it can be seen to form a cycle as described by a closed stress-strain hysteresis loop (Fig. 1a). For the case where successive stress points are either converging or diverging, the hysteresis curves do not form a closed loop (Fig. 1b). For this case the authors assumed that fatigue damage could be attributed to each successive range as half-cycles.
The method was further developed by Okamura et al. [5] and Downing and Socie [6] as a vector based algorithm which identified full RF cycles and half-cycles based on a three-point criteria without the need to rearrange the data series, and enabled efficient utilisation in computer software. This greatly reduced the data http storage requirements as the stress signal could be read into the algorithm in real-time and processed directly into RF cycle spectra. This definition of the algorithm has been refined and included in the ASTM cycle counting standard [2]. Amzallag et al. [7] conducted a wide ranging industry consultation and defined a standardised algorithm which identified RF cycles based on a four-point criterion. The three and four point versions of the algorithm were shown to identify the same cycles by McInnes and Meehan [8], who presented a series of fundamental properties of RF counting to demonstrate the equivalence of the two methods. Although various forms of the RF algorithm exist, the four-point algorithm presents the most unambiguous criterion for the identification of closed hysteresis loops, and is defined below.

Four-point Rainflow counting criterion
RF counting requires the time history to be first processed into a Peak-Valley (PV) series consisting of local maxima and minima which define the turning points, or load reversals, of a time series. Point x m is identified as a local maxima or minima within a time series of length M if, Once the data have been filtered according to the PV criteria, full RF cycles are identified in the range formed by points x n to x n+1 if they meet the four-point criterion, jx nÀ1 À x n j P jx n À x nþ1 j 6 jx nþ1 À x nþ2 j n ¼ 2; 3; 4; . . . ; N À 2 where N signifies the length of the PV filtered series. If the range formed by points x n to x n+1 meets the four-point criterion then the points are recorded before deleting them from the PV series, thus enabling further ranges to be formed between the adjacent points x nÀ1 and x n+2 . The process is repeated until all ranges which meet the four-point criterion are recorded and deleted from the PV series.
Storage of the counted ranges is achieved with a two dimensional histogram to record the cycle stresses. The form of the histogram may be chosen to preserve detailed cycle hysteresis information which may be significant in further statistical analysis, for example with the min-max or max-min matrices where cycles are binned according to the loading sequence [9]. As a minimum, the histogram should record the cycle range and mean stress levels as inputs to final damage calculations.

Rainflow residue
Once all full RF cycles which meet the four-point criterion have been identified and deleted from the PV series, a 'residue' of data points will typically remain. The residue consists of a series of diverging data points from the start to the maximum and minimum points, followed by a converging section of points to the end of the PV data series. Referring to Fig. 2, no remaining closed hysteresis cycles can be identified within a diverging or converging series as no further ranges are bounded by adjacent points of higher and lower value. However, as the stress path formed by the residue constitutes some of the largest ranges in the original series, they should be accounted for if a conservative estimate of fatigue damage is to be made. Two dominant methods exist in the literature to process the RF residue and are outlined in Sections 2.1 and 2.2.
Whenever a subset of a longer time history is RF counted, cycle ranges which are formed between points which span beyond the subset have the potential to be cropped. If there is a large variation in the mean stress level, which is not fully contained within the subset period, then some of the largest cycles will not be accounted for. These cycles are termed ''transition cycles" or ''ground cycles" [10], and a degree of artificiality will be introduced if the residue data points are processed as an isolated set, as closed hysteresis cycles cannot be formed. The only way to accurately identify all RF cycles within a data set according to the four-point criterion is to process the entire time history consecutively, yet the application of RF counting algorithms must always utilise a finite length of data, as chosen by the analyst and by limitations on computational capacity.
Glinka and Kam [11] presented an approach which allowed extended time periods to be read and processed incrementally, thus limiting the required computational capacity by minimising the amount of data required to be handled by the RF algorithm at any one time. A more versatile method is included in Amzallag et al. [7, pp. 292-293] which addresses the same issue by concatenating consecutive residue periods which remain after RF processing. However, although the method allows transition cycles to be accounted for accurately according the four-point criterion, it has not found widespread acknowledgement and no generalised proof of the methodology has been presented.
The three methods of processing the RF residue periods are presented in Section 2 below. An analytical proof is presented in Section 3 which demonstrates the equivalence of cycles which are identified from the residue concatenation methodology outlined in [7] with those which would be identified by RF processing a continuous series. In Section 4, the different approaches are  applied to two case studies in order to illustrate the potential impact of the choice of residue processing method on calculated fatigue damage. Finally, Section 5 outlines the suitability of each of the three methods for practical applications.

RF residue processing methodologies
The three distinct methods available for processing the residue data points are described below and presented in the process diagram in Fig. 3.

Half-cycle counting methodology
This approach is identified in the original definition of RF counting given by Matsuishi and Endo [4], where the authors assumed that each successive range will attribute half a cycle of fatigue damage in the material. From Fig. 2, subsequent half-cycle ranges are identified between points A-B, B-G, G-J, J-O, O-R, R-S, S-V, V-W. At least twice as many ranges will be identified from the residue data points as would be identified as fully closed cycles. Therefore, when the counted residue cycles are stored in the RF histogram the number of cycles added to each bin is reduced by a factor of 0.5. The ASTM RF counting definition of the threepoint algorithm [2, pp. 5-6] is capable of identifying half-cycles which occur up to the maximum data point in the series; after completion of the algorithm, the residue data points following the maximum still remain and must be accounted for as half cycles. Half-cycle counting may be applied directly to the residue which remains from application of the four-point RF criterion, and the resulting cycles can be shown to be identical to those produced by the three-point algorithm.

Simple Rainflow counting methodology
If the stress time history is representative of a repeated loading sequence then all residue data points will ultimately form fully closed cycles as they will fall between repeated extremes. With the four-point algorithm this can be achieved by joining two repeated residues and then reapplying the four-point criterion (Eq. (2)). Closed cycles can then be identified between the repeated maximums, leaving the residue points outside of the maximums which can then be discarded. This is expressed as [7].
From Fig. 2, the residue series is repeated to give a sequence (1) is then reapplied and the repeated point A must be deleted to ensure that the PV sequence is maintained. Eq. (2) is then reapplied to identify closed cycles from all points that fall between J and repeated point O; ranges are formed by points V-W, R-S, B-G, J-O. The remaining points account for the repeated residue, and are therefore discarded.
The simple RF counting methodology is implemented in the three-point algorithm by rearranging the stress time series to start and end with the maximum data point prior to PV processing and Fig. 3. RF counting process diagram for long time periods (modified from [12]). Grey boxes identify steps which relate to the residue processing methods outlined in Sections 2.1, 2.2, 2.3 below. RF counting, and will identify identical cycles [8]. Therefore, the approaches implemented in [2, pp. 6-7], [6, p. 32], [7] and [8] are equivalent.

Residue concatenation methodology
The following steps apply the residue concatenation procedure outlined in [7, pp. 292-293] to the simple case of two PV periods, with reference to Figs. 4 and 5; 1. Define two series of PV processed data points A 1 , B 1 , C 1 , . . ., H 1 and A 2 , B 2 , C 2 , . . . , H 2 . 2. Apply the four-point criterion (Eq. (2)) to both series to identify all full RF cycles. Full cycles are identified between points D 1 -E 1 , and E 2 -F 2 (Fig. 4a). 3. Store the cycles and delete the identified data points D 1 , E 1 and E 2 , F 2 from the respective PV series (Fig. 4b). No more full RF cycles can be identified according to Eq. (2). The remaining points form the two RF residues. 4. Concatenate the two residues in their original chronological order. Apply the PV criteria to the concatenated points H 1 and A 2 to ensure the PV series is maintained. Delete point H 1 (Fig. 4c). 5. Repeatedly apply the four-point criterion to the concatenated series until all fully closed RF cycles have been stored and removed. (Fig. 5a-c). Full RF cycles are identified between points A 2 -B 2 , F 1 -G 1 , C 2 -D 2 sequentially. 6. The remaining residue points are shown in Fig. 5d, and must be processed by either half-cycle or simple RF counting. In practice, successive residue periods may be concatenated to allow additional closed hysteresis cycles to be unlocked.

Equivalence of concatenated series
The residue concatenation method outlined in Section 2.3 can be shown to account accurately for transition cycles between consecutive periods by following the fundamental properties of the four-point criterion presented by McInnes and Meehan [8].
3.1. Identical RF cycles are identified using residue concatenation as would be identified from a continuous series McInnes and Meehan defined an 'end-point bounded sequence' (EPBS) as any series of data points in which the maximum and minimum values lie at the start and end of the series. As all points are bounded between these local extremes, all the ranges contained between the end points will form fully closed cycles according to the four-point criterion. The authors show that the full cycles contained within an EPBS are independent and unaffected by the  data points which lie outside of the sequence (Property 3.5, [8, p. 552]). Furthermore, as the residue remaining from the four-point algorithm consists of a set of EPBSs from the PV series (Property 3.8, [8, p. 554]), all the ranges which meet the four-point criterion will do so independently of the data points which occur before or after the PV series. Therefore, according to these two properties, the specific RF cycles which can be identified from separate periods must also form valid cycles within a continuous series.
McInnes and Meehan also showed that the RF cycles which meet the four-point criterion will do so regardless of the order in which the criterion is applied (Property 3.3, [8, p. 551]). Therefore, the EPBSs which may be formed between residue periods when they are concatenated can release additional RF cycles which would otherwise have been identified from a continuous series.
Referring to Fig. 4a, the data point subsets C 1 , D 1 , E 1 , F 1 and D 2 , E 2 , F 2 , G 2 both form EPBSs and therefore, according to Property 3.5 [8, p. 552], the full cycles formed by the ranges D 1 to E 1 and E 2 to F 2 are independent and unaffected by the concatenation of the two residues. Points C 1 to G 2 then produce an EPBS within the concatenated residues which allows additional full cycles to be formed by the ranges A 2 -B 2 , F 1 -G 1 , and C 2 -D 2 (Fig. 5). According to Property 3.3 [8, p. 551], the same cycles would be identified from a continuous series, although they would have been extracted in a different order.

Deleted residue end points do not affect the correct identification of RF cycles
RF cycles which fall within an EPBS that includes the end point in a residue series are not affected if the point must be deleted after concatenation in order to maintain the PV sequence. This is true because, although the first or last point in an EPBS might be deleted, the sequence would always be extended to a bounding point which is of greater range. Fig. 6 demonstrates the three basic configurations for the connected ends of two concatenated residues (additional permutations exist which are essentially mirror images in the horizontal and vertical planes). Fig. 6a shows the configuration where the PV sequence is maintained and no points are required to be deleted. Fig. 6b shows the case where one point must be deleted to fulfil the PV criteria and the EPBS formed by B 1 to C 1 is extended to B 1 to A 2 as point C 1 is removed. Fig. 6c shows the case where both the concatenated end points must be deleted and the EPBSs formed by B 1 to C 1 and A 2 to B 2 are both extended to the range formed by B 1 to B 2 . Therefore, Property 3.5 [8, p. 552] holds true when RF residues are concatenated and the PV criteria is reapplied.

Residue processing comparison using experimental data
The significance of the different residue processing methods can be compared by investigating their impact on calculated fatigue damage.

Fatigue damage comparison methods
The Palmgren-Miner linear damage hypothesis assumes that the fatigue damage in a loaded component can be expressed as the sum of damages contributed by each stress cycle, where D is fatigue damage fraction, and n i /N i is the ratio of operational cycles to the maximum allowable number of cycles at each  stress range. Although fatigue crack propagation can be found to be influenced by the load sequence, the linear damage sum is commonly used in cases which require a statistical representation of the loading, such as found with environmental loading where the load sequence is rarely well defined at the design stage. Therefore, the linear damage sum has been used for this analysis, although a more detailed fracture mechanics approach may also benefit from the ability to correctly identify large amplitude hysteresis cycles. The maximum allowable number of cycles N is taken from empirical S-N data, as generalised by the Basquin relation, where Dr is either stress range or amplitude (stress range will be used here in accordance with the RF cycle definition), m is the Wöhler exponent, and log a is the intercept of the curve on the log N axis. Combining Eqs. (3) and (4), the ratio of fatigue damages can be expressed as, where n A and n B are the cycle spectra produced by different RF counting methods. It can be seen that, by taking the damage ratio, the log a intercept term cancels and the impact of the different counting methods is affected only by the Wöhler exponent from the S-N curve. It can also be seen that a hypothetical linear increase in stress ranges such as may arise from a stress concentration factor, for example, would also cancel with the damage ratio, indicating that the difference between the RF counting methods would be affected by the underlying load process, but not by the stress magnitude. As a fatigue endurance limit could be exceeded by such a linear increase in stresses any results calculated in this analysis would be trivial; i.e. the use of a constant gradient S-N curve means that the form of Eq. (5) enables the general case to be examined.
S-N data is typically derived from component testing with a zero mean cycle stress. However, stress cycles in the tensile range may produce greater levels of fatigue damage, and S-N test results are known to be strongly dependent on the mean stress level. Mean stress correction models may therefore be used to adjust the stress range prior to damage calculation from S-N curves using, where Dr r is the stress range or amplitude at non-zero mean stress r, Dr 0 is the equivalent stress range or amplitude at zero mean stress, and r UTS is the ultimate tensile strength of the material.
Values of Z = 1 and Z = 2 give the Goodman and Gerber relations, respectively, while the Soderberg relation is given with Z = 1 and r UTS replaced by the yield stress [1]. Mean stress correction models are applicable in the tensile range only, therefore Dr r ¼ Dr 0 is used in compression. The impacts of the residue processing methodologies were then assessed against the following variables; The underlying load processes. The Wöhler exponent used to calculate fatigue damage. The length of the subset period used for the half-cycle and simple RF methods outlined in Sections 2.1 and 2.2.

Measured datasets
Two datasets representing different load processes were used to investigate the residue processing methodologies. The offshore measurement buoy shown in Fig. 7a is governed by dynamic wave loading, but the mooring system is also effected by a semi-diurnal tidal cycle. The offshore wind turbine shown in Fig. 7b also experiences hydrodynamic loading, but the structural response is dominated by the aerodynamic and operational loads.
For the offshore measurement buoy, a load cell was used to measure forces located at the connection between the buoy and a catenary mooring line, recorded with a sample rate of 20 Hz. The force time series was then converted to stresses using the geometry of the attachment shackle. Details of the mooring system are outlined in [13].
For the offshore wind turbine support structure, strain gauge measurements were also recorded at 20 Hz sample rate and converted to time series stresses using the Elastic Modulus of the base material. The gauge was located away from any stress raising features at the base of the wind turbine near the mean water level, axially aligned with the cylindrical tower to measure the nominal bending stresses. Material constants for the measured components are shown in Table 1.

Results
Fig . 8a shows a six month stress history from the offshore measurement buoy. The accumulated fatigue damages were calculated according to Eq. (5), expressed as the ratio of damages produced by RF counting the time series in subsets (using the half-cycle and simple residue processing methods outlined in Sections 2.1 and 2.2) to the damages produced by RF counting a continuous period (using the residue concatenation method outlined in Section 2.3). Fig. 8b shows the damage ratios produced using ten minute subsets; a typical simulation length for a wind turbine load case due to the level of statistical stationarity of the wind loading process found over this time scale [14]. Fig. 8c shows the damage ratios that result using three hour subsets, a length of time which is typically used to characterise wave loading due to assumed sea-state stationarity over this period [15]. Wöhler exponents of m = 3 and m = 5 were used, relating to typical S-N curves for steel components [16]. The data were processed using a verified in-house RF counting code which was developed in MATLAB [17]. Fig. 9 displays a one year period of stress data collected from the multi-megawatt offshore wind turbine support structure and the corresponding accumulated fatigue damage ratios, calculated in the same way. The shapes of the RF cycle spectrums produced by both datasets are shown in Fig. 10.
The final values of the damage fraction ratios, together with the impact of the mean stress compensation models (Eq. (6)), are shown in Table 2.

Discussion
The results presented in Section 4 show that the significance of transition cycles on calculated fatigue damage is dependent upon the underlying load process, the slope of the S-N curve, and the subset length chosen to RF process the data. The catenary mooring system shown in Fig. 8a is representative of the typical scenario where the maximum and minimum stresses in the load process occur well within the time frame which can be processed by conventional RF methods. Longer term variations in mean stress levels do occur over the 12 h tidal cycle, but are of small amplitude in comparison to the stress response from dynamic wave loading. The effect of the varying mean stress is noticeable at the start of the dataset when the transition cycles are more significant in comparison to the low dynamic wave loading. Larger loading events at 2.5 days and 61 days result in closer correlation, with convergence of the calculated damage levels within 97% for each of the RF methods. The longer subset length of 3 h produces a closer correlation between the methods, as more of the tidal cycle is included in the subset period.  The wind turbine data shown in Fig. 9, however, includes a large change in the tower mean stress level as a result of the quasi-static rotor thrust loading which changes with variations in wind speed and direction. The dynamic structural response, which is overlaid on the varying mean stress, is of comparatively low amplitude. The result is that RF counting of the data in subsets using the methodologies outlined in Sections 2.1 and 2.2 accounts for only 37-62% of the fatigue damage which results when the data is processed as a continuous series, using a Wöhler exponent of m = 5. The difference in damage levels are seen to converge after approximately 2-6 months, which is indicative of the length scale of the stress cycles which arise from changes in the quasi-static wind loading. With a Wöhler exponent of m = 3, however, the difference in fatigue damage produced by each of the residue processing methods is negligible, due to the fact that a shallower S-N curve gradient will attribute less weight to the high amplitude/low frequency stress ranges which characterise transition cycles. It should be noted that use of an S-N curve which incorporates a 'knee', or a minimum fatigue endurance limit (which may justify the use of a filter incorporated into the RF analysis to remove stress cycles which are small enough not to effect fatigue life), would increase the impact of the transition cycles because greater significance would be attributed to the high amplitude stress ranges.
The inclusion of mean stress compensation models in the damage fraction ratio results presented in Table 2 is seen to have minimal impact on the significance of transition cycles, with between 2% and 4% less impact with the Goodman and Soderberg relations for the offshore wind turbine dataset with a Wöhler exponent of m = 5. Although the impact is relatively low (the damage from the wind turbine maximum tensile stress range of approximately 30 MPa would be impacted by the Soderberg relation by a factor of (1 À 30/335) À5 = 1.6), this is understandable as cycles from both RF residue processing methods would be effected, and the majority of cycles in the dataset are in compression. The effect of the mean stress compensation models is less than 0.5% in all other cases.
The fatigue damage levels calculated by the half-cycle and simple RF residue processing methods are of comparable magnitudes due to the fact that the bounding maximum cycle range is consistent. Both methods are applicable when the range of stresses produced by the load process is included within the RF counting subset. A typical example is the identification of load cycles experienced by a power station boiler where the load sequence, including maximum and minimum stresses, can be related to long term operational states. The EN standard for Water-tube boilers and auxiliary installations [18] specifies a range of methods for the processing of RF residue points which are comparable to the half-cycle counting and simple processing methods outlined in Sections 2.1 and 2.2.
Concatenation of residue periods, as outlined in Section 2.3, is capable of addressing two main shortfalls of conventional RF methods. Firstly, large volumes of data, which may prohibit RF processing of a continuous series due to computational limitations, can be dealt with correctly. This may be particularly applicable to long term load measurement programmes which typically generate large quantities of data. Secondly, transition cycles which span RF counted periods can be correctly accounted for according to the four-point criterion. Although a factor of two scatter can be expected in S-N test results [19], the identification of RF cycles in a long time history can, and therefore should, be performed accurately.
Additionally, the residue processing method outlined in Section 2.3 could be incorporated into fatigue design by the load case approach through the utilisation of information regarding the long term load history. The design of support structures for multimegawatt wind turbines, for instance, is facilitated by the time domain simulation of structural dynamics under a discretised set of environmental and operational loading conditions [20]. An approach to account for transition cycles with operational wind turbine measurement programmes is outlined in [21], whereby a synthetic time history of maximum and minimum stresses from each load case is constructed corresponding to a long term history of wind speed and direction measurements. However, the approach can be found to be overly conservative as it involves the double counting of data points as both half-cycles within the ten minute load case period and as successive full RF ranges within the synthesised long term stress series. Residue concatenation can avoid this double counting as it accounts for transition cycles correctly according to the four-point criterion.

Conclusions
The three main variations of RF counting described above are distinguished by the way in which the residue is accounted for, and the most suitable method should be selected according to the application. Specific findings include; The concatenation of successive RF residue periods has been shown to enable the same closed hysteresis cycles to be identified as would be produced by RF counting the data as a continuous series. The conventional methods of half-cycle and simple RF counting the residue periods are suitable when the entire stress range seen by a component is contained within the analysed period of data. Concatenation of residue periods is suitable when the data to be processed contains a slowly varying mean stress which results in transition cycles which would otherwise be cropped by half-cycle or simple RF processing the data in subsets. Using this method, transition cycles can be accounted for correctly according to the four-point criterion. The impact of transition cycles is most significant with the use of higher Wöhler exponents.