Fatigue lifetime calculation of wind turbine blade bearings considering blade-dependent load distribution

Rotating bearings are some of the most commonly employed machine elements. As such, they are well-understood and thoroughly researched pieces of technology. Fatigue lifetime calculation is internationally standardized through ISO 281, which is based on the assumption that loads act on a bearing under constant rotation. Blade bearings of wind turbines do not conform to this assumption since their movement typically consists of small, repetitive oscillations. Moreover, their load distribution differs considerably over the bearing circumference, a load case for which ISO 281 refers to ISO 16281 and which requires detailed simulations of the bearing to be sufficiently precise. Aside from ISO 16281, the NREL DG03, a guideline for pitch and yaw bearing lifetime, lists two methods for incorporating bearing loads into the fatigue life calculation. This paper compares all three methods. Two of the methods can not be used directly for the double-row four-point bearing used in this paper and are thus slightly adjusted. Load distributions in the bearing are simulated and curve-fit by means of a novel approach using regression analysis. The method from NREL DG03, which requires the least computational effort, is shown to result in a much higher lifetime than the other two, which are based on internal load distributions of the bearing. The two latter methods are shown to produce very similar results. An adjustment is proposed for increasing the accuracy of that lifetime calculation method which requires the least computational effort in order to resemble the other two more closely.

1 Introduction

Blade bearings are a critical component of any modern wind turbine. Enabling the turbine to pitch can reduce the loads on a multitude of its components significantly. This allows components to have a lighter design and a higher return on investment of the entire machine.

Apart from continuous pitch control (CPC), which turns all blades simultaneously by the same angle, individual pitch control (IPC) has been the subject of comprehensive research (e.g., Bossanyi, 2003; Selvam et al., 2009; Shan et al., 2013). IPC turns blades individually in order to reduce asymmetrical rotor loads, which contribute significantly to fatigue loading of a wind turbine's rotating components (Bossanyi, 2005). However, blade bearings do not necessarily benefit from this load reduction since it is achieved by increased movement of the bearing. Movements with IPC are typically small, repetitive oscillations. This movement pattern differs greatly from that of bearings in most other industrial applications, where bearings usually rotate continuously or, in some cases, turn very seldom at all. Lifetime research has hence mostly been focused on the former, which make up the vast majority of bearings sold and used. Blade bearings thus represent somewhat uncharted territory, which does not sit well with the fact that their replacement is a very costly procedure, which is thus to be avoided.

Bearings exhibit a vast number of possible failure mechanisms, including fatigue, fretting corrosion, brinelling, false brinelling, and more (Stammler et al., 2019). Of these, rolling contact fatigue of the bearing raceways used to be particularly common for rotating bearings. It has therefore been and continues to be the subject of much research (Sadeghi et al., 2009). A breakthrough was achieved by Lundberg and Palmgren (1947), who published a general calculation method for the calculation of rolling contact fatigue. Later, and heavily based on the works of Lundberg and Palmgren, ISO 281 (ISO, 2010b) was published as an international standard for the calculation of raceway fatigue.

ISO 281 is intended for bearings under continuous rotation subjected to a constant axial and radial load or a combination thereof. To account for more complicated load situations, ISO 16281 (ISO, 2010a) was added later, allowing the calculation of fatigue lifetime for any arbitrary load situation.

Even with these standards in place, a great deal of uncertainty remained with regards to the calculation of pitch bearing lifetime. While ISO 16281 allows for consideration of the complicated load situation caused by a tilting moment, oscillatory movement patterns have yet to be considered in any of the standards. Moreover, large slewing bearings behave somewhat differently from the smaller ones on which the standard is primarily based. In 2009, the NREL (National Renewable Energy Laboratory, USA) (Harris et al., 2009) thus published the DG03, a guideline for the calculation of yaw and pitch rolling-bearing life. It collates the state of the art for the lifetime calculation of pitch and yaw bearings and thereby allows for the consideration of the aforementioned factors. However, none of the approaches therein has been verified for large-scale slewing bearings.

Moreover, failure modes of blade bearings are manifold and not just limited to fatigue (see Stammler et al., 2019). Wear is a typical damage mode that commonly occurs due to the small, repetitive oscillations of blade bearings, but no sufficiently reliable calculation methods exist for wear prediction in blade bearings. Hence, life calculation is typically limited to fatigue life. Assessing which failure mode is most common is difficult since, to the authors' knowledge, there are no large data sets of blade bearing failures publicly available.

The present uncertainty is reflected by the certification demands of manufacturers. In its 2003 Guideline for the Certification of Wind Turbines (Germanischer Lloyd, 2004), GL (Germanischer Lloyd) required a rating life calculation for blade and yaw bearings “if applicable”. In its 2010 Guideline (Germanischer Lloyd, 2010), this requirement was removed when the guideline stated that lifetime calculation for pitch bearings was not required for any turbine. Subsequently, in the 2016 Guideline (Germanischer Lloyd, 2016), the requirement was once again changed to require a lifetime calculation according to NREL DG03 under all circumstances.

The present paper examines different lifetime calculation methods from the abovementioned standards and guidelines and compares them to each other in order to highlight differences in the methods and their results. First, the simulations underlying the present paper and the calculation methods used herein are explained in detail. Then, results of the methods presented are compared and discussed.

2 Simulation

The calculation methods presented herein are intended for a double-row four-point contact ball bearing in a nearshore wind turbine. Turbine loads are simulated according to International Electrotechnical Commission (IEC) 61400-1 (IEC, 2019). The load distribution in one of the pitch bearings is simulated by means of a finite-element (FE) simulation that includes the connected blade and hub of the turbine.

2.1 Turbine model

Simulations of the time series were carried out using the IWES Wind Turbine IWT-7.5, a wind turbine model designed by Fraunhofer IWES and described by Popko et al. (2018). It is a nearshore turbine with 7.5 MW rated power output, designed for wind class IEC A1 (IEC, 2019). Additional properties are displayed in Table 1.

Table 1 Main turbine properties of IWT-7.5 (Popko et al., 2018).

Figure 1 FE model of blade, bearing, and hub.

The model assumes that the turbine operates with a controller designed by the German wind turbine manufacturer Enercon. Enercon used the aeroelastic model of the IWT-7.5, equipped with their own IPC controller, to run load simulations. This controller is a wideband IPC, designed to minimize loads as much as possible without limiting the movements of the pitch bearing. The Enercon IPC activates at wind speeds slightly below rated speed. This speed region contributes a large share to the overall fatigue loads of the turbine. The control values are the loads in a nonrotating hub coordinate system, and the main objective is to minimize loads on the steel structures of the turbine (hub, machine frame, tower). Stammler et al. (2019) presented results of different load simulations with another IPC controller for the same turbine. Note that, in contrast to the present work, the controller used by Stammler et al. (2019) only activated IPC above the rated speed.

2.2 Load calculations

Aeroelastic simulations of the wind turbine were carried out according to IEC 61400-1, DLC1.1. A total of 20 years of lifetime were simulated. The calculated loads and movements of one bearing are used as input for the subsequent calculations.

Available literature on the analysis of time series with oscillatory movement patterns is sparse. Fatigue lifetime calculations are most commonly done using a rain flow count (Matsuishi and Endo, 1968), which has proven to be an effective method (see Dowling, 1971). The NREL DG03 also employs a “rainbow cycle” [sic] count in one of its examples. Consequently, this paper uses a rain flow count according to ASTM International E1049 (ASTM International, 2017). Note that for analyzing other types of surface-induced raceway damage, a range pair counting may prove necessary instead (Stammler et al., 2018b). While the sum of all cycles obtained using either a rain flow count or a range pair count will be identical, the length of the cycles will differ. This will have an impact on the calculated lifetime if additional factors for oscillation as used by Schwack et al. (2016) are employed.

After a rain flow cycle count of the bearing oscillations, the results were further divided into bins of the resulting tilting moment, its angle, and the absolute pitch angle of the bearing using the procedure described by Stammler et al. (2018b).

2.3 FE model

The generation of the entire FE model as well as all simulations was performed using Ansys R3. For all FE simulations, a one-third rotor star FE model was used. It consists of a rotor blade, one-third of a rotor hub, a pitch bearing, and a stiffener plate. Using only a one-third rotor star model greatly reduces the computational effort. Doing so makes the model behave symmetrically, meaning that it is not possible to simulate different loads acting on the three blades, a process which is assumed to have negligible effects on the loads of one bearing. Part of the model is shown in Fig. 1. The entire model consists of approximately 855 000 elements and 956 000 nodes, which results in a computational time in the range of 25 min using a PC with an Intel Xeon E5 3.7 GHz processor and 128 GB RAM.

The outer ring of the bearing connects directly to the blade flange of the hub. The bearing’s inner ring connects to the blade with the stiffener plate in-between. Using stiffener plates is a common way to reduce the ovalization of the bearing, which is caused by the blade. The stiffener plate is made of steel and has a thickness of 25 mm, which is a typical thickness for a blade flange of that size. The blade model contains a fully modeled root section. The component contacts are simulated to be bonded over the connecting surfaces, meaning that no bolts and friction-based contact behavior are implemented. Not implementing bolts might lead to a little more flexibility of the bearing rings, which can result in a larger tilting of the rings towards each other. In turn the loads on both rows are distributed slightly less evenly. In-house investigations have shown that this effect only has a very small influence on the bearing's load distribution. Thus, the effect of this simplification is to be assumed negligible. The model is completely fixed at the hub's rotor flange (downwind) and partially fixed at the hub's upwind flange. In addition, cyclic constraints at the hub's one-third cutting planes are implemented. All these boundary conditions enable a realistic deformation behavior of the rotor hub to be modeled. The loads are applied to the blade's spar caps at a blade length of 40 m. Concentrating all loads acting on a blade at one point is a common method for determining bearing loads that does not completely reflect the load application on a real blade. In-house investigations beforehand have shown that doing so delivers comparable results to a realistic load application as long as the load application point is not located in the first quarter of the blade. Figure 2 exemplarily shows the characteristic deformation of the bearing rings which is caused by the surrounding structures and characterized by an imbalance of maximum deformation between the traction and compression side.

Figure 2 Bearing ring deformation in meters for the load case M_y=24 MNm, F_x=600 kN, F_z=1.2 MN.

2.4 Pitch bearing

The bearing used for all simulations is a double-row four-point contact ball bearing as described in Stammler et al. (2019). Table 2 lists the main properties of the bearing. It has the typical design and dimensions of a pitch bearing for a wind turbine of that size. The pitch diameter d_m refers to the distance between two opposite rolling-element centers.

Table 2 Main properties of blade bearing investigated (cf. Stammler et al., 2019).

Figure 3 Bearing cross section and raceway definitions.

The balls of the bearing are not fully modeled but represented by nonlinear springs, which is a common approach according to Daidié et al. (2008). The nonlinear force deflection curve to represent the Hertzian contact between the ball and raceway is calculated according to Houpert (2000). Figure 3 shows a cross-sectional view of the bearing model. The bearing model allows the load distribution between all balls on each of the four raceways to be analyzed.

2.5 Contact forces

During the aeroelastic simulations performed for the lifetime simulations of the turbine, a wide variety of different loads act on the pitch bearings. Of these, the three most influential factors are the resulting tilting moment M, that is, the resulting moment from edge- and flap-wise moments; the angle of said resulting moment, hereinafter called the load angle, β; and finally the pitch angle θ of the blade. These three factors are defined as depicted in Fig. 4. They have the strongest influence on the load distribution in the bearing (Stammler et al., 2018a). While the effect of M as the load onto the bearing is obvious, the influences of β and θ are less apparent. Large slewing bearings such as those used for blade pitching tend to be very elastic. This causes the hub and blade as well as their nonrotationally symmetric stiffness behaviors to change the load distribution in the bearing significantly, depending on their respective orientation. Different load distributions in the bearing then, in turn, create a different equivalent load for the lifetime calculation. The internal axial reaction forces resulting from the tilting moment tend to have a much higher impact on the load distribution within the bearing than the axial force F_z. Nonetheless, a representative value of F_z has been considered for the simulations, depending on the currently acting moment M. Radial forces F_x and F_y are applied at 40 m blade length as described in Sect. 2.3 and thus determined such that they cause the desired tilting moment M at the blade root.

Figure 4 Coordinate system and angle definitions.

Contact forces are thus simulated at a discrete number of points using 358 different combinations of M, β, and θ. Simulations have been run for a grid of data points, shown in Fig. 5. The grid in this case was chosen such that all operating points during the aeroelastic simulations lie within it, hence allowing a regression analysis at all times. Note that, in general, choosing a larger choice of operating points will result in more robust regression analysis results.

Figure 5 All combinations of pitch angle θ, moment M, and load angle β occurring during aeroelastic simulations and the combinations chosen for FE simulations.

2.6 Lifetime calculation methods

All the standards and guidelines mentioned in Sect. 1 calculate rolling contact fatigue lifetime as

$$\begin{array}{}\text{(1)}& L_{\mathrm{10}}={\left({\displaystyle \frac{C_{\mathrm{a}}}{P_{\mathrm{a}}}}\right)}^{\mathrm{3}},\end{array}$$

where L₁₀ denotes the time until 10 % of bearings will show first signs of fatigue on their raceways (Harris, 2001). C_a is the (axial) load rating of the bearing, which for all methods shown is calculated as

$$\begin{array}{}\text{(2)}& C_{\mathrm{a}}=\mathrm{3.647}{b}_{\mathrm{m}}{f}_{\mathrm{c}}{\left(i\cos \mathit{\alpha }\right)}^{\mathrm{0.7}}{Z}^{\mathrm{2}/\mathrm{3}}{D}^{\mathrm{1.4}}\tan \mathit{\alpha },\end{array}$$

according to the ISO 281. C_a is based on geometrical and manufacturing properties such as the row count, i; the number of balls, Z; the ball diameter, D; and the (initial) contact angle, α. The factor b_m equals 1.3, whereas the factor f_c depends on the material and geometry of the bearing and has been chosen according to DIN SPEC 1281-1 (DIN, 2010) for all calculations in this paper. As all these parameters are known at the manufacturing stage, the load rating is generally provided by the bearing manufacturer. P_a of Eq. (1) refers to the dynamic equivalent load, which is a measure for the loads acting on the bearing during operation. Three different calculation methods for determining the dynamic equivalent load are presented in the following, two of which, here referred to as NREL 1 and NREL 2, are listed in NREL DG03.

2.6.1 Method NREL 1

Firstly, an equation based on the applied tilting moment is given as

$$\begin{array}{}\text{(3)}& P_{\mathrm{a}}=\mathrm{0.75}{F}_{\mathrm{r}}+F_{\mathrm{a}}+{\displaystyle \frac{\mathrm{2}M}{d_{\mathrm{m}}}},\end{array}$$

where F_r and F_a refer to the applied radial and axial forces, respectively. M stands for the applied tilting moment and d_m for the diameter of the bearing. Equation (3) does not require any knowledge of the bearing other than its diameter and is thus often used for simple, rough calculations.

2.6.2 Method NREL 2

As a more sophisticated approach, the DG03 also lists

$$\begin{array}{}\text{(4)}& P_{\mathrm{a}}={\left({\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{NREL}}}}\sum _{j=\mathrm{1}}^Z{Q}_j^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}}\cdot Z_{\mathrm{NREL}}\sin \mathit{\alpha },\end{array}$$

which is based on the individual rolling-element loads, Q_j, taken from a simulation of the entire bearing. Strangely enough, Q_j is not further defined despite the fact that pitch bearings tend to be four-point bearings, meaning that up to four forces can act onto a ball at the same time. Equation (4) can thus not be used as it is shown because Q_j could refer to any combination of the four contact forces. This paper therefore uses a slightly modified approach, given by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{P}_{\mathrm{a}}& ={\left({\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{NREL}}}}\sum _{j=\mathrm{1}}^{Z_{\mathrm{NREL}}}(Q_{j\mathrm{A}}\sin \mathit{\alpha }+Q_{j\mathrm{B}}\sin \mathit{\alpha }{)}^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}}\cdot Z_{\mathrm{NREL}}\\ & ={\left({\displaystyle \frac{\mathrm{1}}{Z_{\mathrm{NREL}}}}\sum _{j=\mathrm{1}}^{Z_{\mathrm{NREL}}}(Q_{j\mathrm{A}}+Q_{j\mathrm{B}}{)}^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}}\cdot Z_{\mathrm{NREL}}\sin \mathit{\alpha },\end{array}\end{array}$$

with the loads Q_jA and Q_jB as defined in Fig. 3, where effectively only the axial component of the rolling-element loads (Q_jA sin α+Q_jB sin α) is considered in the calculation. Variable Z_NREL is interchangeably referred to as “the total number of balls in the bearing” and the “number of rolling elements in a row” by NREL DG03, even though pitch bearings are commonly double-row bearings. This paper assumes the former, meaning that $Z_{\mathrm{NREL}}=Z\cdot i$ since otherwise no method to calculate the entire equivalent load for a double-row bearing would be defined.

Note that the variable α is referred to as the “nominal contact angle” in ISO 281 but simply as the “contact angle” without further specification in NREL DG03. For reasons of consistency, all calculations in this paper use the initial contact angle shown in Table 2. Contact angles of a highly loaded pitch bearing will generally increase significantly over the circumference, which would increase the equivalent load as per Eq. (5). However, consideration of a changing contact angle in the calculation of C_a of Eq. (2) (and in the original calculation of P_a according to Eq. 4) is not possible without changes to the overall calculation procedure and has thus not been done in the present work.

2.6.3 ISO 16281

As a third approach, the individual lifetime of each raceway is calculated according to ISO 16281. Like NREL 2, this approach is based on individual rolling-element loads. Equivalent loads Q_ei and Q_ee for a representative contact of the inner and outer ring, respectively, are calculated as

$$\begin{array}{}\text{(6)}& Q_{\mathrm{ei}}={\left({\displaystyle \frac{\mathrm{1}}{Z}}\sum _{j=\mathrm{1}}^Z{Q}_j^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}},\mathrm{and}{Q}_{\mathrm{ee}}={\left({\displaystyle \frac{\mathrm{1}}{Z}}\sum _{j=\mathrm{1}}^Z{Q}_j^{\mathrm{10}/\mathrm{3}}\right)}^{\mathrm{3}/\mathrm{10}}.\end{array}$$

Note that Q_j is now clearly defined as the normal force between one ball–raceway contact. All equations presented in ISO 16281 are intended for single-row bearings, while pitch bearings have two rows with four contacts on each row. Therefore the load ratings Q_ci and Q_ce of the inner and outer ring are calculated as

$$\begin{array}{}\text{(7)}& Q_{\mathrm{ci},\mathrm{e}}={\displaystyle \frac{C_{\mathrm{a}}(i=\mathrm{1})}{Z\sin \mathit{\alpha }}}{\left[\mathrm{1}+{\left\{{\left[{\displaystyle \frac{\mathrm{1}-\mathit{\gamma }}{\mathrm{1}+\mathit{\gamma }}}\right]}^{\mathrm{1.72}}{\left[{\displaystyle \frac{r_{\mathrm{i}}}{r_{\mathrm{e}}}}\left({\displaystyle \frac{\mathrm{2}{r}_{\mathrm{e}}-D_{\mathrm{w}}}{\mathrm{2}{r}_{\mathrm{i}}-D_{\mathrm{w}}}}\right)\right]}^{\mathrm{0.41}}\right\}}^{\pm \mathrm{10}/\mathrm{3}}\right]}^{\mathrm{3}/\mathrm{10}}.\end{array}$$

Refer to ISO 16281 for details on the variables of Eq. (7), and note that C_a(i=1) has been used as ISO 16281 is intended for single-row bearings. These load ratings are then used to calculate the lifetime L_10r,p of each of the four inner–outer raceway contact pairs, p=1…4, according to

$$\begin{array}{}\text{(8)}& L_{\mathrm{10}r}={\left[{\left({\displaystyle \frac{Q_{\mathrm{ci}}}{Q_{\mathrm{ei}}}}\right)}^{-\mathrm{10}/\mathrm{3}}+{\left({\displaystyle \frac{Q_{\mathrm{ce}}}{Q_{\mathrm{ee}}}}\right)}^{-\mathrm{10}/\mathrm{3}}\right]}^{-\mathrm{9}/\mathrm{10}}.\end{array}$$

The raceway lifetime L_10r has thus been calculated strictly according to ISO 16281, aside from the choice of C_a(i=1), which is not explicitly demanded in ISO 16281 but can be reasoned as stated above. Finally, an adjustment of ISO 16281 is undertaken to account for the fact that the bearing has four contact pairs rather than one by calculating the total lifetime of the bearing as

$$\begin{array}{}\text{(9)}& L_{\mathrm{10}}={\left(\sum _{p=\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{10}r,p}^{-e}\right)}^{-\mathrm{1}/e},\end{array}$$

where $e=\mathrm{10}/\mathrm{9}$ denotes the Weibull modulus for ball bearings (cf. DIN SPEC 1281-1; DIN, 2010). Comparing this approach to Eq. (4), it can be seen that the NREL DG03 approach is essentially an abbreviated version of the ISO 16281 using some simplifying assumptions. For informative purposes, the lifetime calculated as per Eq. (9) is then turned into an equivalent load for the entire bearing according to Eq. (1) by

$$\begin{array}{}\text{(10)}& P_{\mathrm{a}}={\displaystyle \frac{C_{\mathrm{a}}}{L_{\mathrm{10}}^{\mathrm{1}/\mathrm{3}}}},\end{array}$$

where $C_{\mathrm{a}}=C_{\mathrm{a}}(i=\mathrm{2})$ now, as usual, denotes the load rating of the entire bearing. This allows the three methods to be compared on the basis of the equivalent loads they provide.

2.7 Additional factors for lifetime calculations

The basic rating life L₁₀ as calculated according to the equations in Sect. 2.6 does not consider a variety of other factors, most notably here the movement patterns of the bearing and its lubrication conditions.

2.7.1 Oscillation

A number of approaches exist for factoring in the movement patterns of bearings (see Schwack et al., 2016, for a comparison). Oscillation affects the fatigue lifetime of a bearing in a number of ways, and particularly the inner (oscillating) ring will be differently loaded than on a continuously rotating bearing. Moreover, other damage mechanisms such as wear are significantly accelerated by small oscillations (see Stammler et al., 2019). With the oscillation angle θ_osc of a rain flow cycle bin measured in degrees, this paper assumes

$$\begin{array}{}\text{(11)}& n=n_{\mathrm{osc}}\cdot {\displaystyle \frac{{\mathit{\theta }}_{\mathrm{osc}}}{\mathrm{180}}}\end{array}$$

for all calculation methods, meaning that the various effects of oscillation on fatigue lifetime are not considered. This approach is mathematically equivalent to that presented by Harris (2001) (cf. Schwack et al., 2016).

2.7.2 Lubrication

Low rotational speeds worsen lubrication conditions, which, in turn, reduces the fatigue lifetime of a bearing. ISO 281 proposes an approach based on the multiplication of the basic rating life L₁₀ with a factor a_ISO to receive the modified rating life

$$\begin{array}{}\text{(12)}& L_{\mathrm{10}m}=a_{\mathrm{ISO}}\cdot L_{\mathrm{10}},\end{array}$$

which is also recommended by NREL DG03. Some differences exist in the calculation of a_ISO, but mostly the approach given by NREL is adapted from the ISO under the assumption of poor lubrication conditions. Effectively, for the simulations given, these differences have very little influence on the results. The temperature has been assumed to be 15 ^∘C, and the corresponding viscosity of the lubricant was thus calculated according to DIN 51563 (DIN, 2011) based on two points at 40 and 100 ^∘C. The value did not exceed a_ISO=0.1003 at any operating point or with any method, while its minimal possible value is defined as a_ISO=0.1 since a calculation for lower values is not possible according to the current state of knowledge as per ISO 281. Exchanging the ISO approach for that of the NREL thus only shortened the lifetime by about 0.15 % for simulations done in this paper.

3 Results and discussion

First, the approach for contact force regression analysis used herein is verified. Thereafter, the results of all three lifetime calculation methods listed in Sect. 2.6 are compared. An adjustment for the simplest method of all three is given so that it more closely resembles the other two.

3.1 Contact force regression analysis

As depicted in Fig. 5, there are significantly more data points of the aeroelastic simulation than there are FE-simulated points. A novel approach is presented to determine the contact forces Q_j for all points of the aeroelastic simulation. A regression analysis of the form

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{Q}_j& =\left[c_{M,\mathrm{0}}+\sum _{k=\mathrm{1}}^{k_{\max }}{c}_{M,k}{M}^k\right]\cdot \\ & \left[c_{\mathit{\beta },\mathrm{0}}+\sum _{l=\mathrm{1}}^{l_{\max }}\left(c_{\mathit{\beta },s,l}\sin \left(l\mathit{\beta }\right)+c_{\mathit{\beta },c,l}\cos \left(l\mathit{\beta }\right)\right)\right]\cdot \\ & \left[c_{\mathit{\theta },\mathrm{0}}+\sum _{m=\mathrm{1}}^{m_{\max }}\left(c_{\mathit{\theta },s,m}\sin \left(m\mathit{\theta }\right)+c_{\mathit{\theta },c,m}\cos \left(m\mathit{\theta }\right)\right)\right]\end{array}\end{array}$$

is used for each ball–raceway contact j in the bearing. Variable k_max refers to the degree of the polynomial used to approximate moment M, while l_max and m_max denote the weights for the Fourier series used to approximate pitch and load angles θ and β. Once Eq. (13) is expanded, it is linear with regards to the various combinations of its c variables. These combinations can then be determined by means of a least-square fit using results from the FE simulations. The problem being a linear least-squares problem, its solution can be determined with low computational effort using, for instance, the Moore–Penrose inverse of the problem. Even though most of its summands are not necessary for the regression analysis, the problem can thus be used as shown because the computational time remains short.

Equation (13) essentially consists of one factor for the moment, load angle, and pitch angle, respectively. While the moment is approximated with a polynomial, both angles are approximated with a Fourier series. To consider interdependencies of the three factors, they are multiplied by each other.

Figure 6 shows interpolated forces of one contact j for a fixed pitch angle of 10^∘. Changes in the load angle β locally take on a sinusoidal shape. For certain angles, the force disappears completely as the balls lose contact. An increase in the moment M unsurprisingly leads to higher contact forces.

Figure 6 Forces Q_j for a contact j at 0^∘ for a fixed pitch angle of $\mathit{\theta }={\mathrm{10}}^{\circ }$.

Likewise, the pitch angle affects the contact force as shown in Fig. 7 for another contact j and at a moment of M=10 MNm. As the blade pitches, stiffness properties of the inner bearing ring relative to the hub-fixed coordinate system change. This causes a different load distribution, which, for the contact shown, results in lower forces for a higher pitch angle. Contacts at other positions will behave differently.

Figure 7 Forces Q_j for a contact j at 0^∘ for a fixed moment of M=10 MNm.

The degree of M and the weights of θ and β should be chosen with the number of FE data points in mind. Note that with a grid as shown in Fig. 5, particular care has to be taken since in areas where the resolution is less dense (as for high moments in this example) the approximation may behave differently than in those where the resolution is higher (as for low moments in this example). Choosing k_max=3, l_max=2, and m_max=2, the average error per contact in the simulated positions is then about 1100 N.

After regression analyses have been done for each Q_j, a contact force distribution is determined for each operating point as shown in Fig. 8 for the example of M=20 MNm, $\mathit{\beta }={\mathrm{90}}^{\circ }$, and $\mathit{\theta }={\mathrm{10}}^{\circ }$. Results from the FE simulation are shown with solid lines; those from the regression analyses performed are plotted with dashed lines. Since the calculations performed in Eqs. (4) and (6) are a type of weighted average over all contact forces, small differences in the FE simulations have negligible effects.

Figure 8 FE-simulated (solid lines) and interpolated (dashed lines) forces Q_j for M=20 MNm, $\mathit{\beta }=\mathrm{90}{}^{\circ }$, and $\mathit{\theta }=\mathrm{10}{}^{\circ }$.

3.1.1 Lifetime calculation methods

A bin counting is carried out in order to calculate the lifetime of the bearing. For 54 different bins of the oscillation angle, the three variables moment M, load angle β, and pitch angle θ are put into 24, 70, and 90 bins, respectively. In total, each oscillation angle bin is thus, in turn, divided into $\mathrm{24}\times \mathrm{70}\times \mathrm{90}=\mathrm{151}\mathrm{200}$ bins. For each of these bins, the frequency of occurrence is calculated separately. The result is shown in Fig. 9, where point sizes reflect the frequency of occurrence. Not surprisingly, most of the bearing operation takes place at low pitch angles, with $\mathit{\theta }=\mathrm{0}{}^{\circ }$ being the most common one. Moreover, load angles between 0 and $-\mathrm{180}{}^{\circ }$ are also rare as these represent wind coming from behind the turbine. Consequently, it pitches out of the wind at these operating points. The maximum moment of M=24 MNm is only achieved at a load angle close to $\mathit{\beta }=\mathrm{90}{}^{\circ }$ and a pitch angle close to $\mathit{\theta }=\mathrm{10}{}^{\circ }$.

Figure 9 Dots for each chosen class occurring during aeroelastic simulations, with size relative to frequency of occurrence.

For each of the bins, the equivalent load is determined according to the three variants laid out in Sect. 2.6. Moment M unsurprisingly has the strongest influence on the load. Figure 10 shows the calculated equivalent load for all existing bins with a pitch angle of $\mathit{\theta }=\mathrm{10}{}^{\circ }$. The overall ratio between P_a and M can be seen to be almost linear, as assumed by variant NREL 1 (see Eq. 3). However, the calculated loads differ between the three variants. The NREL 2 and ISO 16281 variants are strikingly similar, with the NREL 2 being slightly higher in magnitude. This can be attributed to the fact that the NREL 2 method is a simplified, slightly more conservative version of the ISO 16281. Compared to the other two, variant NREL 1 is much lower: at the highest loads for M=24 MNm, this results in a difference of 12 % compared to the other two. As the equivalent load is factored into the lifetime to the power of 3 (see Eq. 1), this difference will have a considerable effect on the calculated lifetime, especially considering the fact that it occurs near a common operating point of the turbine as shown in Fig. 9.

Figure 10 Equivalent loads P_a for a fixed pitch angle of $\mathit{\theta }=\mathrm{10}{}^{\circ }$ at operating points from the aeroelastic simulations.

With a constant moment of M=5 MNm, the effects of the load and pitch angle can be seen in Fig. 11. As already observed above, the NREL 1 variant produces lower results than the other two. At a value of $P_{\mathrm{a}}=\mathrm{2.48}\times {\mathrm{10}}^{\mathrm{6}}$ N, it is at least 13 % lower than the other two variants are at their lowest points. This difference will have a strong impact on the calculated lifetime. Moreover, it is constant for all values as it does not consider the effect of changes in θ and β (see Eq. 3). The other two variants, however, are based on internal load distributions of the bearing. Changes in θ and β are thus reflected by a change in equivalent load P_a. As can be seen, the reaction to these differences occurs with different sensitivity. For all existing bins shown, variant NREL 2 results in the highest equivalent load P_a since assumptions made during its derivation from the ISO 16281 are mostly conservative in nature. The difference remains below 10⁵ N at all times, which equals about 3 % of the maximum load for this load case. The qualitative behavior is similar as well: for both variants at $\mathit{\theta }=\mathrm{25}{}^{\circ }$, a minimum is reached at $\mathit{\beta }\approx -\mathrm{20}{}^{\circ }$. This corresponds to the situation where a moment is acting edgewise on the blade, which is why the spar caps are not carrying any load while the maximum pressure is acting on the side of the hub, which has softer stiffness behavior. With an increase in the load angle β at the same pitch angle $\mathit{\theta }=\mathrm{25}{}^{\circ }$, a maximum is then reached at $\mathit{\beta }\approx \mathrm{100}{}^{\circ }$. The moment is now primarily acting flap-wise, and the spar caps are carrying most of the load, which causes the blade to exhibit stiffer behavior. Furthermore, the spar caps are pushing into the downwind side of the hub, which similarly exhibits stiffer behavior. The overlap of the high blade stiffness due to the spar caps carrying most of the load and the stiffer backside of the hub resisting against this pressure cause the highest resulting load to occur in this position. With an increase in pitch angle θ, one can see that the stiffness behavior of the blade significantly impacts the equivalent load P_a as the position of the maxima and minima changes with θ. At a pitch angle of $\mathit{\theta }=\mathrm{70}{}^{\circ }$, the maximum is then closer to load angle $\mathit{\beta }=\mathrm{70}{}^{\circ }$ since that will be the direction of a flap-wise moment.

The difference between the maximum and minimum of the two methods is then approximately 5×10⁵ N, which is roughly 14 % of the maximum load shown. This difference reflects the impact of θ and β. It will impact the calculated lifetime to some extent but not nearly as significantly as the resulting moment M.

Figure 11 Equivalent loads P_a for a fixed moment of M=5 MNm at operating points from the aeroelastic simulations.

Differences in equivalent loads between NREL 2 and ISO 16281 decrease with an increase in the moment, as shown in Fig. 12 for M=20 MNm near the highest foreseen moment of the bearing. Note that, once again, only bins that actually occurred during aeroelastic simulations are shown. Maximum differences in equivalent loads between NREL 2 and ISO 16281 now reduce to 0.25 %, and differences caused by the load and pitch angles are at most in the range of 7 % since the turbine does not change them significantly when the resulting moments are high. This range will be the most significant for the lifetime calculation since it occurs frequently and does so under high loads.

Figure 12 Equivalent loads P_a for a fixed moment of M=20 MNm at operating points from the aeroelastic simulations.

Other operating points not shown in Figs. 10 and 11 remain similar with respect to their qualitative behavior. The NREL 1 variant remains lower than the other two for every single bin examined, and equivalent loads of variant NREL 2 are higher than those of ISO 16281 in 99.9 % of cases.

Using these values, the overall lifetime of the bearing L_10m is calculated, paying due consideration to the frequency of occurrence of each bin. The results are shown in Fig. 13. Additional factors have been chosen as explained in Sect. 2.7, where the lubrication parameters have been used according to ISO as differences to the methods presented in NREL DG03 are negligible (see Sect. 2.7.2). As expected, the considerable differences of NREL 1 compared to the two other methods examined here have a strong impact on the calculated lifetime. Using variant NREL 1 for the calculation of equivalent loads hence leads to roughly 1.7 times the lifetime of the other two methods since the calculated loads are lower. Fatigue lifetime is therefore predicted to be 4273 h, or roughly 178 days. The other two methods barely differ with regards to their results. Both predict a lifetime of around 107 d. This resemblance can be attributed to the similarity in their equivalent loads at common operating points of the turbine, as seen in the figures above.

These calculated lifetimes are significantly below the expected lifetime of a wind turbine of 20 years. This discrepancy is similar to other publications on the issue (cf. Schwack et al., 2016), and even calculated examples in the NREL DG03 remain well below 20 years after consideration of all factors. While available data on blade bearing failures are sparse, the calculated lifetimes are so low that no data are necessary to disprove them; the mere fact that blade bearing exchanges are costly, time-intensive operations suffices to illustrate that they will typically have to last longer.

Assessing the exact reasons for the differences between calculation and reality is, however, difficult. Calculation methods are largely based on research with small bearings, whose conditions during manufacturing and operation differ from those of large slewing bearings (cf. Göncz et al., 2010) such as pitch bearings. Large slewing bearings will, in relation to their size, generally deform more than small bearings. The stiffness of the connected structures such as the blade and hub will highly affect the deformation behavior and lead to large changes in the contact angle of 20^∘ and more in a highly loaded state (cf. Chen and Wen, 2012), which is not considered in any calculation approach used in this paper. The usage of a_ISO≈0.1 may also be questioned since it reduces the calculated lifetime to a tenth of its initial value. Moreover, the lifetime L_10m denotes the statistical point in time at which first damage occurs on the raceway for 10 % of bearings. This view of a lifetime might be too conservative for pitch bearings, which have to be as slender as possible to enable a high return on investment of the entire wind turbine and thus may continue to be operated when they are already damaged. While small bearings are expected to fail very soon after the first surfaced raceway damage, large bearings may be more robust in this regard. Lastly, as mentioned in Sect. 2.7.1, various effects of the oscillatory movement of pitch bearings are not considered in the results of this paper.

Figure 13 Fatigue lifetimes of all methods investigated.

3.1.2 Adjustment of NREL 1

Given the fact that the qualitative behavior of NREL 1 closely resembles that of the other two methods, as can be seen in Fig. 10, and the fact that pitch angle θ and load angle β influence the lifetime less significantly than the resulting moment M, Eq. (3) provides a good basis for a simplified lifetime calculation. In order for NREL 1 to result in a similar lifetime, specifically the term $\frac{\mathrm{2}M}{d_{\mathrm{m}}}$ should be adjusted as it represents the strongest influence on the resulting equivalent load P_a. For the simulations in this paper, an adjustment to

$$\begin{array}{}\text{(14)}& P_{\mathrm{a}}=\mathrm{0.75}{F}_{\mathrm{r}}+F_{\mathrm{a}}+{\displaystyle \frac{\mathrm{2.5}M}{d_{\mathrm{m}}}}\end{array}$$

generates lifetimes as depicted on the right-hand side of Fig. 13. The result of adjusting NREL 1 is thus slightly lower than that of its two counterparts, thereby allowing for some margin of error stemming from changes in pitch and load angle if different operating points were to occur during simulations. The calculation is, however, much simpler and thus well suited for rough analyses of the raceway fatigue lifetime of blade bearings. The results of all methods are compared in Table 3. Equation (14) is thus valid for the specific turbine examined in this paper. The authors assume that for a different turbine or, more specifically, a different combination of hub, bearing, and blade, a different adjustment may be appropriate.

Table 3 Fatigue lifetimes of all methods investigated.

4 Conclusions

Blade bearings of wind turbines operate under unusual operating conditions compared to others in the industry. Some details which go beyond the internationally standardized calculation of fatigue lifetime as per ISO 281, such as the calculation of equivalent loads or the consideration of oscillatory behavior, can thus be obtained by a number of methods. This paper investigated three different approaches for the calculation of equivalent loads P_a required for the lifetime calculation: two according to NREL DG03 and one according to ISO 16281.

For the case of a blade bearing of the reference wind turbine, load distributions in the bearing have been simulated and curve fit to allow for consideration of a variety of operating points. The results show that changes in the load and pitch angle of a rotor blade bearing lead to significant changes in the equivalent load P_a. However, the impact of the resulting moment was identified to be more significant than that of the load and pitch angle.

The two methods that calculate P_a on the basis of simulated load distributions (NREL 2 and ISO 16281) have been shown to provide very similar results. The third method (NREL 1), which is merely based on global loads acting on the bearing and which does not require detailed simulations of the latter, has been shown to result in much higher lifetimes than the other two methods. An adjustment has been proposed for NREL 1 to match the results of its two counterparts more closely.

As already seen in other publications on the fatigue lifetime calculation of pitch bearings, these results are far lower than the expected lifetime of a turbine of 20 years. Many reasons may account for this discrepancy: calculation methods are largely based on research with small bearings, whose conditions during manufacturing and operation differ from large slewing bearings. Effects such as the changing contact angle during operation that occurs for large slewing bearings with flexible attached structures are not considered at all. Moreover, the effect of a_ISO, which reduces the calculated lifetime to a tenth of its initial value, may be put into question. On top of that, the definition of the lifetime L_10m, which denotes the statistical point in time at which first damage occurs on the raceway for 10 % of bearings, may be too conservative for blade bearings.