Abstract
Fatigue life calculation in bearings is an issue that has been the focus of decades of research culminating in national and international standards that are easily applicable to a wide variety of applications. However, almost all research has focused on simple, repetitive operating conditions. This paper explains why these approaches are not completely applicable to bearings under certain conditions, particularly not to those with a change in the load position over time undergoing inconstant movement. Based on existing standards, namely ISO 281, the paper then proposes a numerical approach to calculate the life of a bearing under stochastic operating conditions. The bearing life for any arbitrary load history and movement patterns can be calculated using the presented approach. The approach is presented and some example applications are given, notably including the life calculation of a rotor blade bearing in a wind turbine.
1 Introduction
Fatigue life calculation in bearings is a much studied topic that has been the focus of over a century of research [1]. Various different models have since been proposed, see Refs. [2,3] for an overview. Most influentially, Lundberg and Palmgren [4] published a model that continues to form the basis for most, if not all, national and international standards, most importantly ISO 281 [5]. Almost all of the research deals with bearings subjected to repetitive operating conditions. Their load is mostly constant or varies merely in magnitude, and their movement tends to be rotational or, in rare cases, oscillating with a constant oscillation angle that keeps repeating itself [6,7].
Bearings in some industrial applications do not adhere to these assumptions. Section 2 explains when these assumptions do not hold true. Notably, some bearings in wind turbines, such as yaw and pitch bearings, experience stochastic loads and movement patterns. Their load varies depending on wind conditions, not just in magnitude but also in direction. The movement of the bearing is determined by the pitch controller of the turbine, which adjusts the blade angle to control the power of the turbine and to reduce fatigue loading on other components of the turbine [8]. This leads to loads and movement as shown in Fig. 1, which is very much unlike anything seen in most typical industrial applications. Leupold et al. [10] proposed a segmentation of the bearing, which could presumably be used for such movement and loads, but use empirical S–N curves derived for merely one particular bearing.
This paper proposes a numerical approach to calculate the fatigue life of such bearings based on ISO 281 [5] and thus also Lundberg and Palmgren [4]. Based on a scalable rolling contact fatigue model from the two aforementioned references, a bearing is segmented. The cumulated damage of each segment over all operating conditions experienced by the bearing during its operating time is calculated first, and only thereafter are the individual segment damage cumulations combined to obtain one bearing life. This approach allows local changes in load magnitude and differences in the number of load cycles to be captured properly.
Some effects that can be modeled with the proposed approach have been considered in other publications. This includes some publications on oscillating bearings, notably Breslau and Schlecht [6] and Houpert and Menck [7], who derive analytical equations for an oscillating bearing in one invariable operating condition. The proposed approach in this paper can be considered as a generalized version of ISO 281 and the aforementioned publications. It hence leads to similar or identical results if identical operating conditions are assumed, but it can also be used for a wider variety of other conditions, most notably arbitrary timeseries of load magnitude, direction, and movement of the bearing.
The approach is explained in Sec. 3. Thereafter, some example applications are given in Sec. 4. These include a very simple load case, which is identical to the assumptions used by ISO 281 [5] and consequently leads to the same result, since the approach is derived from it. It also includes more advanced operating conditions which cannot be calculated purely with ISO 281, which are compared with other available literature. The Appendix gives a derivation of proportionality constant k required for Sec. 3 based on Lundberg and Palmgren [4] and ISO 281 [5].
2 The Palmgren–Miner Hypothesis Applied to the Entire Bearing
The Palmgren–Miner hypothesis was initially suggested by Palmgren in 1924 [11] to combine the effect of varying load levels in bearing life calculations mathematically. It has since found utilization in various applications, far exceeding just the use for bearings. However, the assumption that fatigue damage at all considered load levels occurs in the same position is implicit to the theory. As will be illustrated in the following, the application of the theory to an entire bearing life can therefore lead to erroneous outcomes.
2.1 Derivation of Palmgren–Miner Applied to the Entire Bearing.
2.2 Error/Simplification of Palmgren–Miner Applied to the Entire Bearing.
Equation (1) and, subsequently, Eqs. (2) and (3) implicitly assume that the load distribution as well as the load cycles occur in the same position over all considered load levels. Moreover, they assume that the shape of the load distribution remains constant.
Lundberg and Palmgren [4] use a bottom-up approach in their calculation of fatigue life. In the event of uneven loading over the circumference of the bearing, as found in radial bearings or bearings subjected to a bending moment, they state: “The probability that the whole [···] ring will endure is the product of the probabilities that the different sections of the ring will endure”.1 Consequently, they take account of variable loads first on a local, segment-wise level and then combine individual segment lives into a total one2 in their derivation of L.
This very local consideration of loading is then what Eq. (1) is incapable of taking into account, since it works with lives Li of an entire bearing. Variable L gives the life of a bearing in one constant condition, but can no longer give information about local effects of it. Operational changes that lead to a vastly different load situation can thus not be correctly captured. Figure 2 depicts two load scenarios that lead to the same life L. If the bearing condition were to change, causing it to operate in either scenario 50% of the time, Eqs. (1)–(3) would predict no change as compared to a constant load direction. This is obviously incorrect and it is a consequence of the loss of local information in the derivation of a bearing life L. This particular example has little impact on the life of a small, rotating bearing where the inner ring determines the life of the entire bearing, but for a large oscillating bearing, the influence is stronger.
Likewise, an oscillation that changes its mean position cannot correctly be calculated with the approach shown in Eq. (1). As the mean angle of the oscillation changes, different parts of the raceway are loaded, but the local information about this loading is lost in the calculation of a complete bearing life L. The NREL DG03 [19] recommends to correct for oscillation in its Eqs. (17)–(19), and then applies Palmgren–Miner to the entire bearing in its Eq. (22), even though the mean angle of oscillations will change in a typical blade bearing, thereby leading to the aforementioned problem.
Publications based on the DG03 as Schwack et al. [20] and Menck et al. [9] also apply Palmgren–Miner to the entire bearing, and thus suffer from related shortcomings. Similarly, Wöll et al. [21] propose a “Numerical Approach” for irregular operation conditions that applies Palmgren–Miner to the entire bearing and hence does not consider local effects properly.
Depending on the application, the influence of this effect varies. In many cases, the application of the Palmgren–Miner rule to the entire bearing will still give results very close to a more precise calculation, and in those cases, it can be considered a reasonable simplification. For other cases, notably for oscillating operation conditions or large bearings where the inner ring impacts the total bearing life more, the difference to a more precise calculation as shown in Sec. 3 may be stronger.
3 A Numerical Fatigue Life Model
3.1 Equivalent Loads Per Segment m.
The novelty in the approach described herein lies in the application of the Palmgren–Miner rule for each segment over all operating conditions that occur over the operating time of the bearing. Only thereafter are the individual cumulated segment damages combined into a total bearing life L.
The time series may represent a condition that repeatedly occurs for x times, as is common in wind turbine simulations. In this case, the so-called multipliers x are applied to each load cycle n such that xn = x in a given time series. To this end, Eq. (9) needs to include multipliers in the values of as well as in the sum .
3.2 Combining the Segments.
3.3 The Life L.
3.4 Further Details.
In addition to the above mathematical derivations, some details remain to be defined.
For the determination of load cycles, the model analyzes the ball movement relative to each of the rings. For example, a movement θi of the inner ring with a stationary outer ring leads to a movement θc = 0.5 · θi · (1 − γ) of the cage.5 For a stationary outer ring, this can be seen to be the movement relative to the outer ring; however, relative to the moving inner ring, the movement of the cage is given by θi − θc. For each time-step of the simulation, and depending on the position of the cage, each ball’s location with respect to the inner and outer segmentation is then respectively analyzed. No slippage of the rolling element set is assumed. The starting position of the first ball is chosen to be 0 deg, with the rest of the balls being evenly spaced accordingly.
A load cycle then occurs according to Fig. 4 for the definition used in this paper: it takes place in that step in which a segment is left. It does not matter on which side it entered, and it does not matter how long it resided in the segment, or if it resided at all. Therefore, shear stress τ0 and depth z0 used for some cycle n in Eq. (9) are given by the corresponding load Q in that time-step in which the segment is being left.
This definition of a load cycle does not follow directly from the Lundberg–Palmgren theory and one might also choose it differently, but the author expects differences from another load cycle definition (e.g., when a ball enters a segment) to be negligible for a sufficiently fine segmentation.
As in Lundberg and Palmgren [4] and ISO 281 [5], the order of loading is neglected by the given approach.
Moreover, the empirical exponents c, e, h, and the empirical parameter A used in the derivation of Lundberg and Palmgren [4] have to be defined. This publication follows the choice of ISO [5] which are mostly based on Ref. [4]. The values used herein are shown in Table 1 and the Appendix.
Lastly, Eq. (6) requires the variables τ0, z0, and a. For reasons of consistency, this paper employs the methods also used in Ref. [4]. The Hertzian contact is calculated by iteratively solving elliptical integrals by Hertz as given in Lundberg and Palmgren’s Eq. (101), and the Hertzian contact is then defined by their Eq. (41). The variables z0 and τ0 are then obtained by solving their Eqs. (9) and (10). However, there are various other simplified methods to determine these values in the literature, including one by Houpert [22] and another one by Brewe and Hamrock [12]. These may be used as well, if desired, and will lead to similar, albeit not identical, results as the solution of the respective elliptical integrals. It is important to maintain the same Young’s modulus and Poisson’s ratio for the calculation of contact pressure and variable k in Eq. (A6).
4 Example Applications
Some example calculations are illustrated in the following. The approach can be used for any kind of movement in a bearing subjected to any kind of load history.
4.1 Simple Axial Bearing.
For the basic case of a rotating, axially loaded bearing under a constant load and with a contact angle of 90 deg, the calculation can be simplified extremely. In this case, both the inner and outer rings are loaded identically. For a rotating bearing, the amount of load cycles can also be taken from the literature instead of being simulated as per Fig. 4. The amount is generally given by u = Z/2 · (1 ± γ) load cycles per rotation, where Z is the number of balls and the upper sign (+) refers to the inner ring while the lower one (−) refers to the outer ring. For a contact angle of 90 deg, we obtain γ = 0. Thus, each position on the inner and outer ring experiences the same amount of load cycles u = Z/2.
We will use a ball bearing defined by Z = 147, pitch diameter dm = 4,690 mm, ball diameter Da = 80 mm, and groove conformity f = 0.5319. The bearing is loaded with Fa = 10 MN, which is assumed to spread evenly on the rolling elements,6 giving a contact force of Q = 68.03 kN for each ball. Iteratively solving the Hertzian elliptical integrals and applying further equations by Ref. [4], this results in τ0 = 559.57 MPa and z0 = 0.749 mm, a = 9.45 mm, and . With one rotation as reference, we obtain N = u. The constant k = 1.441 · 10−42 according to Eq. (A8).
Since the entire ring is evenly loaded and experiences the same amount of load cycles in every position, it is sufficient to use one 360 deg-wide segment.7 Inserting the above values into Eq. (6) then returns ln(1/Sm) = 8.672 · 10−7. Due to only using one segment, it follows from Eq. (10) that ln(1/Sr) = ln(1/Sm). For the entire bearing, Eq. (11) gives ln(1/SB) = 2 · ln(1/Sr). From Eq. (13) it follows that ξ = 20,192. Since the input data used corresponds to one revolution, it follows that L10 = 20,192 rev.
Applying Eq. (A6) for the above bearing, we obtain Ca = 2.723 MN. The life calculated according to Ref. [5] thus results in L10 = (Ca/Fa)3 · 106 = 20,192 rev, which is identical to the result calculated above. This is due to the simple load case considered, which corresponds exactly to the assumptions used by Lundberg and Palmgren [4] and, by extension, ISO [5].
4.2 Oscillating Bearing.
Oscillating bearings have been considered analytically in the literature based on modifications of the Lundberg–Palmgren theory. Notably, Breslau and Schlecht [6] and Houpert and Menck [7] published approaches to correct the standard ISO life for oscillation. The same is possible with the numerical model presented here, as it is presented in Sec. 3.
To this end, simulations have been performed to reproduce the results of Houpert and Menck [7] using the model presented in this paper. An entire ball bearing, consisting of an inner and outer ring,8 has been simulated. A 7220 type bearing with 15 balls of diameter Da = 25.4 mm and a contact angle of α = 45 deg was used to this end. Inner and outer ring were each divided into 1,800 segments. One oscillation of θ degrees amplitude was performed using the same angle definition as in the aforementioned reference, shown in Fig. 5, for a variety of θ. At the same time, a load zone as in the given reference was used (identical to that one used by Lundberg and Palmgren [4]) and is described by the value ɛ. Simulations have been performed for the same values of ɛ as in the reference to obtain L10, osc. For each value of ɛ, one rotation was performed as well to obtain L10, rot. Then, an oscillation factor was calculated according to , which was then further divided by the so-called Harris factor aHarris = 90/θ for improved display.
The results are shown in Fig. 6. When comparing Figs. 6–10 of Houpert and Menck [7], the factor given here is below that shown by Houpert and Menck [7] for small oscillations. This is because the simulations performed in here automatically include an effect not yet considered in the given reference: With oscillations below the critical angle θcrit = 2π/(Z (1 ± γ)), only a part of the raceway is loaded. This leads to an additional effect on the life, which Houpert and Menck [7] introduce only after showing their9 Fig. 10. That is to say the simulations performed with the model shown in this paper automatically include in their results all relevant effects that occur when any kind of oscillation, small or large, is performed. The user does not require advanced understanding of any effects that arise for any specific type of movement, neither do they have to derive lengthy equations for any number of effects that arise; as long as the actual movement and loads are included in the input data for the model, their effects are automatically included.
4.3 Application of the Methods to a Rotor Blade Bearing.
Rotor blade bearings are subject to stochastic external conditions from the wind, and their load and movement thus cannot be described analytically. An example of rotor blade movement and loads is given in Fig. 1. This is hence where the approach herein can be used to calculate a correct solution based on the assumption made by Lundberg and Palmgren [4] as shown in Eq. (5). The method described in Sec. 3 is thus applied to DLC1.1 time series data given in Ref. [24] and determined according to IEC 61400-1 [25]. This contains simulated load series based on stochastic wind conditions that are assumed to be Weibull distributed according to IEC [25] for wind turbine class IA. The same time series have also been used in the life calculation in Menck et al. [9]. The data are made up of 50 ms time-steps and each time-step is analyzed for load cycles resulting from movement of the bearing, which are then added according to Sec. 3.1. The bearing, too, is identical to that one in the aforementioned publication. It is a double-row four-point bearing with four inner and outer raceways, respectively. The raceways are called 1A, 1B, 2A, and 2B. For a sketch of the bearing, see Ref. [9].
The life of the bearing is calculated according to Sec. 3.2. The contact forces Q are calculated based on finite element simulations according to Menck et al. [9], and τ0, z0, and a are calculated according to Lundberg and Palmgren [4]. Variable k = 9.1396 · 10−43 according to Appendix A. A convergence analysis was performed and resulted in <0.1 % change for more than M = 1800 segments, which is thus the number used for the simulations.
Twenty years of operation have been simulated with 255 10-min bins under different wind speeds and with multipliers corresponding to the respective percentage of the wind speed in question. As a measure of damage, the value ln(1/Sm) after 20 years of operation is displayed in Fig. 7. Each dot in the figure corresponds to a segment. The results of two segments right next to each other can vary considerably. This is particularly due to the large number of balls in each row of the bearing (Z = 147) and the small movements that it performs during operation, which means balls tend to stay within a few segments.
The resulting life is L10 = 858.54 days. This can be compared to the result obtained in Menck et al. [9], which gives a modified reference life of 2520 h or 105 days. includes the lubrication factor aISO ≈ 0.1, correcting for which results in a life of L10 = 1050 days. This is 22% above the result calculated herein, particularly because Menck et al. [9] did not correct for the influence of oscillation on the life. However, even though the aforementioned publication uses Palmgren–Miner applied to the entire bearing which, as stated in Sec. 2, is not correct, the results are reasonably similar. One influence on the relatively small difference is the Weibull slope of e = 1.1 according to Table 1, which is very close to 1. Other authors [1,23] have suggested values for e that differ more from 1, which would result in higher differences between the life calculated from a “global” application of the Palmgren–Miner hypothesis as shown in Sec. 2 and the more accurate “local” application of the finite segment method as shown in Sec. 3 and as used for the results shown here. Lastly, despite the change in the direction of the resulting moment, one can see in Fig. 7 that most of the damage occurs merely within a range of about 60 deg on the compression and the traction side, which means that the change of load direction is not as big as, for example, depicted in the more extreme example of Fig. 2 for this particular case.
5 Conclusions
The paper showed that the application of the Palmgren–Miner hypothesis to the lives of entire bearings can lead to incorrect results if the qualitative position of the load changes, either by a change of the global load direction or by small, nonrepetitive movement that leads to load cycles occurring in different locations. A numerical model was proposed that accounts for variations in the local load cycle position and magnitude for arbitrary time-varying movement and load angles. For the same assumptions as taken by Lundberg and Palmgren [4] (e.g., constant load direction, full rotations, same load zone, Hertzian rolling elements) and consequently ISO 281 [5], the model has been shown to produce identical results to ISO as well as a publication on oscillating bearings by Houpert and Menck [7]. The model was then applied to a rotor blade bearing which undergoes effectively random small oscillations and load changes due to the wind. The effect of this was analyzed and compared to an earlier publication by Menck et al. [9] that applied Palmgren–Miner to the entire blade bearing. Even though doing so is not entirely accurate, as shown in this publication, the results are reasonably similar.
Footnotes
Also known as Weibull weakest link theory according to Weibull [18].
Note that this is an extremely simplified description. For survival probability S, Lundberg and Palmgren [4] start off with . In the case of a pure thrust load, the raceway is evenly loaded and a segment-wise analysis is not required. Hence, in that case, V refers to the entire loaded volume of the raceway. Merely for uneven load distributions on the standing (typically outer) ring, they take the product of all individual survival probabilities, with V → dV. In rotating bearings, this is not needed for the rotating (typically inner) ring (even in the case of an uneven load distribution), since all segments of it experience the same load eventually.
V is used here as , with l referring to the raceway circumference. Lundberg and Palmgren [4] merely define V ∼ az0l, and this paper uses a proportionality constant 1; this proportionality constant may be chosen arbitrarily as long as the constant k (here derived in the Appendix) is derived using the same definition.
This publication uses parameters c, e, and h as given in DIN ISO [5] and derives k based on the same reference. However, other parameters could equally be used. The same principle could also be used for different life models as long as they are based on the Weibull weakest link theory.
γ is defined in the Nomenclature.
It does not spread evenly in real applications, of course, hence why ISO includes λη in Eq. (A6).
Of course, multiple segments could be used as well and would lead to the same result. Likewise, simulating the actual movement of the bearing and counting load cycles as described in Sec. 3.4 would (approximately) lead to the same result as assuming N = Z/2 for one rotation.
All calculations in this paper include the inner and outer ring. This is only emphasized here because Ref. [7] also includes a separate calculation for each ring.
The effect is included in Houpert and Menck [7], Section 4.4, with the factor in their Eq. (45), which has to be applied to the oscillation factors given in their Fig. 10 in order to obtain Fig. 6 of this paper. Further slight deviations may occur because Houpert and Menck [7] use the model by Dominik [23] and adjust it for point contact, giving very similar but not identical results to the model by Lundberg and Palmgren [4], and because possible interdependence of various effects that occur is more accurately captured in the model presented in this paper.
LP use kg and mm as units, in which case A = 10.
Internal calculations not included in this paper have produced about 10% discrepancy in the calculated life without the term in question. As shown in Sec. 4.1, inclusion of the term virtually removes any difference to the standard for simple reference load cases.
Acknowledgment
Thanks to Iman Talai for implementing the model outlined in an earlier version of this paper, and thanks to Matthias Stammler for proofreading the paper.
Funding Data
The German Federal Ministry for Economic Affairs and Climate Action (Grant reference numbers 0324344A and 03EE2033A).
Nomenclature
- a =
half major contact ellipse axis, mm
- c =
empirical exponent
- e =
Weibull slope
- f =
groove conformity
- h =
empirical exponent
- i =
index of an operating condition
- k =
proportionality constant, for other variables using N and mm
- l =
number of revolutions; also: raceway circumference, mm
- m =
index of a finite segment
- n =
index of a load cycle; also: rotational speed, 1/min
- p =
load-life exponent
- q =
percentage of operating time, %
- r =
index of a raceway
- u =
load cycles per rotation
- x =
time series multiplier
- A =
empirical parameter
- C =
load rating, N
- D =
accumulated movement, deg
- I =
number of operating conditions
- L =
life, (106) revolutions
- M =
number of finite segments
- N =
number of load cycles
- P =
equivalent load, N
- Q =
rolling body load, N
- R =
number of raceways
- S =
survival probability
- T =
duration of time series, s
- V =
loaded volume, mm3
- Z =
number of balls
- aosc, aHarris =
corrective factors for oscillation
- dm =
pitch diameter, mm
- dV =
loaded volume of a finite segment, mm3
- lcoll =
collective number of revolutions
- z0 =
depth of τ0, mm
- Da =
ball diameter, mm
- Dn =
raceway diameter, mm
- Fa =
axial force, N
- Lcoll =
collective life, (106) revolutions
- L10 =
life with S = 0.9, (106) rev, days
- L10m =
modified reference life, (106) rev, days, h
- Mres =
resulting bending moment, MNm
- Pcoll =
collective equivalent load, N
- SB =
bearing survival probability
Greek Symbols
Appendix: Derivation of k
The following derivation is heavily based on the original work of Lundberg and Palmgren [4] (termed “LP” in the following) and cannot be completely understood without it. The equations that follow thus use the same nomenclature as LP and their variables are not explained in the following. Furthermore, most explanations given in the original work are skipped here. The focus lies on using the same equations as LP, but with inclusion of a variable k from the very beginning in order to obtain an equation as shown in Eq. (6) rather than merely a proportionality as given in Eq. (5). Harris and Kotzalas [12] give an overview over the calculation with similar, albeit not identical, nomenclature and may be used to follow the equations as well.
Equations from the original work of Lundberg and Palmgren [4] are labeled “LP.” Differences due to the derivation of k are highlighted in (see color equations online).
A.1 Derivation of Basic k
A.2 Derivation of A1
In the following, the derivation of the load rating C will be undertaken in order to determine A1. Following the derivation by LP, the constant A1 is simply replaced by A = 0.089 · A1, with A = 98.0665 according to DIN [26] if N and mm are used as units.10 However, further additions made at the end of the derivation of the load rating C are not included in the aforementioned definition of A1. Particularly ISO [5] has made several adjustments to the load rating following the original derivation by LP, and in the following we will use A1 to include any and all adjustments made to the load rating up until its final form as given in ISO 281.
A.3 Final k
Conflict of Interest
There are no conflicts of interest.
Data Availability Statement
The authors attest that all data for this study are included in the paper.