Method to determine the local load cycles of a blade bearing using flexible multi-body simulation

Conventional methods for designing rolling bearings against fatigue assume that a bearing ring is fully rotating and that the load is ideally distributed over the rolling elements. Blade bearings in wind turbines, are operated under oscillating motions and dynamic loads. The load distribution is strongly dependent on the stiffness of the bearing rings and the surrounding structural components. This has been shown in numerous studies using FEM simulations for static load cases. In this paper a method is presented that reduces the calculation effort of the deformation of the bearing rings, so that a flexible integration into an aeroelastic mbs model of a wind turbine is possible. Thereby an average accuracy of 6.5% between FEM and mbs could be achieved. The model allows the determination of time series of the global load distribution of each raceway. By data processing of the simulation results, the number of load cycles and the maximum contact pressure for individual segments of the raceways could be determined and their fatigue probability could be estimated using the linear damage hypothesis according to Palmgren-Miner.


Introduction
Blade bearings in wind turbines connect the hub and the rotor blade, so that a rotation of the blade is possible and thus the power input and occurring wind loads can be controlled. Since the blade bearing is also part of the supporting structure, a failure of the blade bearing can lead in the worst case to a loss of the rotor blade.
In the majority of wind turbines of the 3 MW classes, double-row four-point contact bearings are used as blade Fig. 1 Cross section of a two row four point bearing. The load is transferred in the marked diagonals for axial force and bending moments bearings. As shown in Fig. 1, the bearing consists of two rows of rolling elements and eight raceways. Depending on its direction the load is transmitted between the inner and outer ring via four possible load diagonals. The angle at which the loads are transmitted changes with increasing load. Under normal conditions only one diagonal of each rolling elements is loaded.
In comparison to conventional bearings, blade bearings are not operated in full rotation, but stand still for a large part of the operating time or oscillating with small angles in the range of approx. 0.2-10°. Additionally, high dynamic bending moments occur due to wind and weight forces [1]. Fig. 2 Overview of the used tool chain to simulate the local load cycles Due to the unusual movement pattern and the resulting poor lubrication of the rolling contact, most common damages of the rotor blade bearings are related to wear, e.g. false brinelling [2]. However, control strategies such as individual pitch control (IPC), which increase the movement activity of the bearing, can also bring fatigue damages of the raceway to the fore [3].
In conventional approaches for blade bearing design, the cross-section loads and pitch movements are determined by a aeroelastic multi-body simulation (mbs) model of the turbine. Then the occurring oscilltions movements are counted and classified by the oscillation angle [4,5]. The load is averaged during the oscillation movement. The load direction and load cycles due to dynamic loads while the bearing is standing still are usually neglected. The influence of the stiffness of the surrounding structural components such as hub and rotor blade on the global load distribution on the individual rolling elements is well investigated for static load cases but usually only simplified considered during the design of the bearing against fatigue [5][6][7][8][9].
In this paper, a method of determining the global load distribution through flexible mbs and local load counting to determine the load cycles of each raceway segment is presented. Thereby, the change of the global load distribution due to the deformation of the bearing rings, the distribution of the load on different raceways, the dynamic load and load direction are taken into account.

Methods
An overview of the presented method for determining the local load cycles on individual sections of the bearing raceways over the entire service life of the wind turbine is shown in Fig. 2. In a first step, an aero elastic mbs wind turbine model [10] is extended by a flexible hub and bearing rings. Using force elements [11], time series of the global load distribution on each rolling element can be calculated, under consideration of the deformation of the surrounding structural components. All relevant design load cases according to IEC 61400 [12] can be simulated. The time series of the global load distribution are converted to a local contact pressure distribution using Hertzian contact theory. Finally, the pressure cycles occurring on each single segment over time are counted using a load counting algorithm and extrapolated to the total service life of a turbine.

Flexible mbs wind turbine model
The starting point of the calculation of the global load distribution is a conventional aeroelastic model of a turbine, modelled in Simpack. The general design parameters of the wind turbine and the blade bearing are given in Table 1. To determine an accurate global load distribution in the bearing, the flexibility of the bearing rings and the surrounding structural components must be taken into account. For this purpose, two separate parts are modelled in FEM, the hub connected the outer ring and the inner ring connected rotor blade. Additional 27 retained nodes with three translatory degrees of freedom are defined on each bearing row and connected to a local area of the raceway. To be able to use the bodies in the mbs environment, both are modally reduced using the Craig-Bampton method [13]. In contrast to FEM, the deformation of the bodies in mbs is calculated using constraint modes and constrained normal modes. The constraint modes represent the static deformation of the retained nodes under load. They are calculated by fixing all but one degree of freedom of the retained nodes and determining the static shape of the deformation. Thus, the number of constraint modes corresponds to the number of degrees of freedom of the retained nodes. The constrained normal modes represent the natural vibration of the body. They are calculated by fixing all degrees of freedom of the retained nodes and calculating the eigenmodes of the body. In Simpack deformation caused by external forces and constraints are evaluated only for the defined retained nodes. Thereby, the degrees of freedom are reduced enormously and thus the computation time. The number of retained nodes of the reduced model was selected in such a way that a good agreement with simulation results of static load cases in the FE model is guaranteed. In the mbs environment Simpack the individual rolling element loads are calculated using the force element 88: "Rolling Bearing". The stiffness of the rolling contact is calculated using the Hertzian contact theory and modelled with non-linear springs. As shown in Fig. 3, for every time step the displacement of each spring is determined by the two closet flexible nodes of outer and ring. The calculated normal force is applied locally on these flexible nodes. The load is distributed to both nodes on the inner and outer ring using a spline interpolation. This allows larger relative rotation between both bearing rings and thus a realistic representation of the pitch movement without load jumps on the flexible bearing rings.
With the described mbs model, time series of the bearing load distribution and the position of the rolling elements are calculated for the most fatigue relevant design load case DLC 1.1, according IEC 61400. [5]. For each average wind speed in the range of 3-25 m/s six different wind filed with a length of 800 s are simulated.

Analytical contact model
Based on the Hertzian contact theory [14] the width (w) and the contact pressure distribution p(x) of the contact ellipse can be calculated for the simulated normal forces of each rolling contact. An increase of the maximum contact pressure due to a truncation of the contact ellipse at high loads is neglected. The geometry parameters of the bearing are given in Table 1. As shown in Fig. 4, each raceway is divided into 12,371 individual segments leading to an element width of 0.58 mm. Depending on the position of the rolling element, the local pressure is applied on specific areas of the raceway. Both the motion of the inner ring and the rolling elements are taken into account. It can be assumed that the rolling elements migrates over time. This means that at a certain pitch angle, not always the same raceway segments are loaded. In order to take this into account, the pressure distribution calculation is carried out for different starting positions of the rolling elements in the range between 0-2.56°which corresponds to the ball distance.

Local load counting
Based on the time series of the local contact pressure, the number of overrollings is counted for each segment. By means of a bin counting, the number of load cycles are determined and classified according to the maximum contact pressure. The interval width of the bin counting is 25 MPa. Thus, the absolute frequency for different maximum contact pressures can be determined for each segment. Load cycles due to dynamic loads when the bearing is not moving are neglected in this this work, but can be taken into account with the presented method without much additional effort. The amount of load cycles determined for the different wind speeds are extrapolated for a wind turbine lifetime of 20 years based on the probability of occurrence for an IA wind location. The probability of occurrence of the individual wind speeds is determined by means of the Rayleigh function Pr given in DIN 61400. Pr To estimate the fatigue behaviour, the damage equivalent load (DEL) is determined for each segment. The different load collectives are summarized using the Palmgren-Miner linear damage hypothesis [15].
Where i is load case number, j the total number of loads, ni the number of load cycles and Ni number of load cycles to failure.
The number of load cycles to failure Ni, is determined based on a fatigue test carried out on a full-size blade bearing and a slope of a single logarithmic S-N curve of -0.33, which is usual for ball bearings [16]. Since the fatigue be- Fig. 4 Calculation of the local pressure according to the Hertzian contact theory havior varies a lot within a bearing type, a large number of tests are necessary for reliable results. Currently, not enough tests have been carried out to determine a reliable service life. The damage factors presented in this paper are only given to clarify the method.

Flexible mbs wind tubrine model
Comparison of static load cases In a first step simulation results of the global load distribution with the conventional modelling approach in FEM and the simplified modelling in mbs for static load cases are compared. Therefor a model consists of a rotor blade, stiffening ring, blade bearing and hub is built up. The hub is fixed at the flange to the main shaft in all degrees of freedom. At the rotor blade in 20 m distance to the blade bearing forces in flap-and edgewise direction are applied, so that for each load case a resulting bending moment of 8 MNm is applied to blade bearing. These loads correspond to the maximum load during power production determined in the aeroelastic mbs model. The load direction is varied between -45°, 0°and 45°and the pitch angle between 0 and 35°. Fig. 5 shows, the global load distribution of the FEM and MBS model for six different load cases. Additionally, the theoretical load distribution for rigid bearing rings is shown. The colors represent the four possible load diagonals, which are shown in Fig. 1. The influence of flexible bearing rings depends on the load direction and pitch angle. The influence is very high for a load angle of -45°with a maximum difference of 100% between the FEM model and rigid system. The reason therefor is a comparatively low stiffness in the area around of 320°. This results in a high reduction of the load distribution in this area. For a load angle of 0°, the influence is with a difference of 25% much lower, due to a comparatively homogeneous stiffness in the highest loaded area. The load distributions simulated with the reduced mbs model show a very good agreement with the results of the FEM model for all load cases. The average error is 6.5% considering all rolling elements above a load of 10 kN. The maximum positive deviation is 5.6 kN, the maximum negative is -3.8 kN. The difference of the service life time according to ISO TS 16281 [17] is about 1.5% between the FEM and mbs model and about 56% between the FEM and the rigid model.
The simulation time is reduced from about 1.2 h to about 2 sec. Thus a sufficient model reduction could be achieved to integrate the blade bearing into the aeroelastic wind turbine model.

Load series for the turbine model
The mbs wind turbine model is able to simulate time series for the global load distribution for the four load diagonals and the pitch angle.
In only one of the three rotor blades the blade bearing is implemented. The simulation time increases by a factor of 4 compared to an aeroelastic mbs model with flexible rotor blades and tower considering natural frequencies up to 20 Hz. Fig. 6 shows an example of the load distribution of the load diagonal 12, according to Fig. 1, for six different time steps, of a simulation of a wind field with an average wind speed of 14 m/s. Visible is the angular shift of the load maximum due to a changing load direction and the increase of the maximum normal forces due to an increased amplitude of the applied bending moment. In addition, the influence of the stiffness of the structural components on the load distribution is shown by their deviation from an ideal sinusoidal shape.

Local pressure distribution
Using Hertzian contact theory, the simulated normal forces are converted into local pressures and assigned to specific segments of the eight different raceways, based on the pitch angle-dependent position of the rolling element set and bearing rings. Fig. 7a shows the contact pressure for all segments of raceway OR 12 according to Fig. 1, for two different time steps. Due to the dynamic load and direction, the shape of the load distribution varies and different segments of the raceway are loaded due to a different pitch angle. Fig. 7b shows the time series of the contact pressure for two different segments of the same raceway. It can be seen that the level of the contact pressure and the number of load  cycles at the segment at 180°is greatly increased compared to the segment at 60°. The temporal succession of the load cycles is strongly dependent on the pitch activity shown in Fig. 7c. This shows oscillations in the range of 0-20°.

Local load counting
First, the maximum pressure and number of load cycles for each surface segment is determined for each of the 136 time series. The individual time series are weighted based on their probability of occurrence of its individual average wind speeds defined in DIN 61400 and extrapolated for a turbine lifetime of 20 years. Fig. 8 shows the absolute frequency of load cycles as a function of contact pressure for each segment of the two raceways of the upper bearing row of the outer ring. The main wind direction is at 180°. Due to the dominant bending moment caused by the thrust, both raceways are loaded mostly one sided. The size of the loaded area is approx. 205°for raceway OR 11 and 260°for raceway OR 12. The highest contact pressure occur on both raceways with a slight offset for the main wind direction.
In order to estimate the service life of the bearing, the damage factor is determined according to Eq. 2. The damage factors of each individual segment of the raceways is shown in Fig. 9. All of them show a clear highest loaded zone and thus a clear preferred location of fatigue can be identified. This results from the dominant bending moments with a clear preferred load direction. However, the load zones are wider than the load distribution of static load cases shown in Fig. 5, which results from the consideration of the shifting load direction. Similar to the load distributions of the static load cases, it is noticeable that the individual raceways are loaded differently and thus there is a different risk of fatigue.
Naturally, the normal forces for the inner ring and outer ring of one load diagonal are similar. The raceways of the outer ring are concave and those of the inner ring are convex, so that the contact pressures for the same normal forces are slightly higher. Since the diameter of the bearing is much bigger then the ball diameter this effect is comparatively small. Due to the movement of the inner ring up to 30°, the load is distributed over a larger area compared to the outer ring, which leads to a slightly reduced number of load cycles. These two effects lead to a different damage factor distribution of both raceways. Notable is a large overlap of the damage factors of the raceways of one bearing row in the range of 50°and 250°. For assembly reasons, there is a filler plug at one position of the raceway to fill the raceways with rolling elements. This is a weak spot at the raceway and should therefore be as unloaded as possible. This is not possible if the load zones overlap strongly.

Conclusion
In this paper a method was presented to simplify common FEM modeling approaches for the simulation of global load distribution in blade bearings, so that an integration into an aeroelastic mbs wind turbine model is possible. It is able to calculate time series of the global load distribution for all raceways of the bearing. Thereby an average accuracy of 6.5% could be achieved compared to the modeling in FEM.
By further data processing the number of load cycles for individual areas of the blade bearing could be determined. It could be shown that the load distribution function for the individual segments is very different. Additionally, the areas with the highest risk for fatigue damage could be identified. The modeling method could be used in the long run to design rotor blade bearings against fatigue.
Load cycles due to dynamic loads with stationary blade bearings have not been considered so far. This could also be taken into account in subsequent work. In addition, fatigue criteria such as Fatemi-Socie [18] or Dang Van [19] could also be applied on the simulation results.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.