Probability assessment of static overload in wind turbine blade bearings considering turbulence, design, and manufacturing variability

This study presents a probability assessment of static overload of a double-row, four-point contact ball blade bearing under ultimate limit state conditions. The methodology follows structural reliability principles applied to this specific limit state. The National Renewable Energy Laboratory 5 MW reference wind turbine is selected for the study, and a Monte Carlo simulation is used to assess static overload reliability and estimate the probability of failure. The sensitivity of the probability of failure to uncertainties in turbulence intensity, material properties, and bearing dimensions is evaluated. Within the bounds examined, the blade bearing conformity has the largest effect on the probability of failure, and ball diameter is next. The probability of damage for case-study wind sites around the world is assessed, and it is observed that the probability of failure is higher in some cases than for wind conditions described by IEC 61400-1.

1 Introduction

Demand for wind turbines will grow significantly over the next decade, and wind turbines are one of the main keys to meeting the IEA Net Zero Emissions by 2050 Scenario (IEA, 2023). The efficient operation of these turbines depends on various components, with blade bearings, also named pitch bearings, being a critical element. Blade bearings facilitate the smooth rotation of the turbine blades, allowing them to capture the kinetic energy of the wind at different wind speeds without significantly increasing the structural loads on the turbine. Blade bearings serve as the connection point between the blades and the hub, allowing the blades to rotate around their axis. Although the entire pitch system assembly costs less than 1 % of the wind turbine (Stehly et al., 2024), changing a blade bearing is costly due to the need to lower the blade with a large crane. The cost of replacement with a crane can reach up to USD 350 000 per week (Mishnaevsky and Thomsen, 2020).

A common and distinguishing feature of the blade bearings is that they involve a rather slow oscillatory motion (Harris et al., 2009). This movement pattern differs from that of bearings in most other industrial applications, where bearings usually rotate continuously (Menck et al., 2020). Unlike most oscillating bearings, blade bearings perform stochastic oscillations rather than constant amplitudes back and forth (Stammler et al., 2024). They oscillate through only a few degrees, up to 20° (Keller and Guo, 2022) during normal power operation, depending on wind conditions and turbine sizes, and up to 90° in emergency feathering, rather than completing full revolutions. The limited rolling distance, long dwell periods at fixed pitch angles, and frequent load reversals with high axial offsets concentrate stress cycles within a small contact zone on the raceways. Blade bearing failure consists of different damage modes, including rolling contact fatigue, core crushing, edge loading, ring fracture, rotational wear, fretting, false brinelling (Andreasen et al., 2022), and cage damage. Many of these failure mechanisms initiate in high-stress regions associated with maximum ball loads. While this study focuses on static overload, applied loads and bearing parameters that drive this mechanism are also relevant to other failure modes.

As part of the design and certification process, the static safety factor of the blade bearing under the ultimate limit state (ULS) must be assessed, as mentioned in IEC 61400-1 (2019), DNV-ST-0361 (2016), Harris et al. (2009), Germanischer Lloyd (2010), and Stammler et al. (2024).

Several studies have analyzed blade bearings. Among them, Menck et al. (2020) studied different lifetime calculation methods from the standards and guidelines and compared them to each other, highlighting differences in the methods and their results, and Schwack et al. (2016) compared different fatigue lifetime calculation methods and showed huge differences in the calculated fatigue lifetime for different approaches. There are not many studies analyzing the static safety factor in blade bearings. Keller and Guo (2022) studied the static load rating and safety factor of the blade bearing of a 1.5 MW wind turbine. They compared the ISO 76 (2006) methodology with the National Renewable Energy Laboratory (NREL) pitch and yaw bearing design guideline (DG03) (Harris et al., 2009) and concluded that the DG03 methodology is the same as the ISO 76 (2006) recommendation for applications subjected to shock loads. Rezaei et al. (2023) studied the blade bearing of the 5 MW NREL reference wind turbine and assessed the variation in blade bearing fatigue with the shear power law exponent, turbulence intensity, and even each individual turbulent wind time series. In another work (Rezaei and Nejad, 2023), the fatigue life of the blade bearing is compared in different wind sites and with IEC-designed blade bearings. The results show that fatigue life at some wind sites is lower, even though their average wind speeds fall below the IEC 61400-1 (2019) category thresholds.

Even though it was shown that the results of wind sites are not the same as those of IECs, it is not clear how reliable the results are, as wind is a stochastic phenomenon. Haus et al. (2024) from more than 55 GW of wind plant data show that blade bearings installed pre-2016 perform fairly well, only reaching a 10 % replacement rate in 15 years. However, blade bearings installed post-2016 on larger wind turbines are projected to have a 10 % replacement rate in only 7.5 years.

The design of blade bearings is governed by a combination of turbine-level and component-level standards. IEC 61400-1, clause 9.8.4, requires that blade bearings demonstrate a minimum static safety factor against permanent deformation in the ultimate load cases. This requirement is based on limiting the local contact stress between the balls and raceways to a threshold value. DG03 expands on this by recommending a contact stress limit of 4200 MPa, based on ISO 76, and defining a methodology to evaluate loads, contact geometry, and bearing strength. DG03 is widely used in the industry but is expected to be superseded in the coming years by the proposed IEC 61400-18, which will standardize pitch and yaw bearing design. While IEC 61400-8 (2024) provides structural design guidance for nacelle and hub components, it does not explicitly include blade bearings within its scope. Nonetheless, it may become a valuable reference if adopted in future standards such as IEC 61400-18, just as IEC 61400-4 (2025) references it for gearbox structural components.

This paper analyzes the probability of static overload in a blade bearing at the ultimate limit state using a structural reliability framework, with a deeper focus on the effect of the wind conditions – particularly turbulence intensity – as well as uncertainties in bearing loads, material strength, and manufacturing tolerances. The goal is to quantify the probability that static contact stress exceeds a specified limit, not to assess all failure modes. While static overload is used here as a representative limit state due to its stress threshold (4200 MPa), the methodology is generalizable. Extreme loading, combined with uncertainty in geometry and material properties, is relevant to other critical failure modes – such as ring cracking – if the appropriate limit state definition is available. A sensitivity analysis is also performed to identify which parameters (e.g., raceway conformity, ball diameter) most strongly influence the probability of static overload. Bearing static damage is considered a criterion for the ultimate limit state. ISO 76 (2006) states that experience shows that a total permanent deformation of 0.0001 of the ball diameter at the most heavily loaded contact does not impair the subsequent operation. In the present study, we treat contact stresses that reach this deformation limit as attaining the ultimate limit state; i.e., stresses approaching the ISO 76 threshold are assumed to represent the onset of static failure and therefore increase the probability of bearing failure. This permanent deformation can cause stress concentrations of considerable magnitude and the formation of cavities in the raceways. These indentations, together with conditions of marginal lubrication, can also lead to surface-initiated fatigue damage (Harris and Kotzalas, 2006). Although some research has shown that the slewing bearing can tolerate higher total permanent deformation while no core crushing occurs (Stammler et al., 2024), those results might be correct for a certain range of bearings with specific material and heat treatment (Lai et al., 2009). The current work assumes a contact stress of 4200 MPa as the criterion for static overload; however, there is a possibility that indentation and core crushing damage do not occur in all bearings at this level.

2 Case study: wind turbine and wind sites

2.1 Reference wind turbine

Load calculations were carried out using the NREL 5 MW reference wind turbine by Jonkman et al. (2009). It is an offshore wind turbine, designed for wind class IEC IB. The turbine properties are displayed in Table 1.

Table 1 The 5 MW NREL reference wind turbine specification (Jonkman et al., 2009).

2.2 Blade bearing

The blade bearing model is considered from work by Rezaei et al. (2023). The bearing is a double-row, four-point contact ball bearing with a total of 250 balls. The general specification of the bearing is presented in Table 2. The details of the bearing specifications and dimensions can be accessed at Rezaei et al. (2023).

Table 2 Blade bearing main dimensions.

2.3 Wind sites

The wind regimes consist of IEC category wind fields and wind sites. IEC category wind fields consist of three wind speed classes: I, II, and III. Each wind speed class has four subclasses, including A+, A, B, and C, based on turbulence intensity. The basic parameters of the wind turbine classes are presented in Table 3 (IEC 61400-1, 2019). In the table, V_ave is the annual wind speed, V_ref is the reference wind speed averaged over 10 min, and I_ref is the reference value of the turbulence intensity.

Table 3 Basic parameters for an IEC category wind turbine (IEC 61400-1, 2019).

The wind sites include 13 in Iran (SATBA, 2022), Pakistan (World Bank Group, 2023b), Vietnam (GIZ, 2023), Ethiopia (World Bank Group, 2023a), Denmark (Ørsted, 2022), and the United States of America (Jager and Andreas, 1996). The wind data include the mean and standard deviation of 10 min wind speed. In order to account for the wind's seasonal effect, data need to cover an entire year. The nominated wind sites cover a whole year's measurements. More information on the wind sites is presented in the work of Rezaei and Nejad (2023).

In this study, the extreme turbulence model (ETM) was investigated. ETM is calculated according to IEC 61400-1 (2019) and prescribes rarer, higher-turbulence realizations that are used later in our ULS analysis. The extreme turbulence intensity at wind sites refers to the maximum turbulence intensity within each wind speed range. To connect the IEC references to the site data, the ETM turbulence intensity curve for IEC classes was computed and compared with the per-bin extreme TI observed at representative sites. Figure 1 shows that Aysha and Kebribeyah exhibit extreme TI values substantially above ETM for class IA (and, at some wind speeds, even above IA+), while Thanh Hai, Mil Nader, and Flatirons are closer to ETM.

Figure 1 Comparison of turbulence intensity for IEC ETM and extreme conditions at selected wind sites.

The extreme turbulence intensities at the wind sites are presented in Table 4. The effects of high turbulence intensity in some wind sites, such as Sujawal, Aysha, and Kebribeyah, are shown and discussed in the Results section.

Table 4 Extreme turbulence intensity in wind sites.

3 Methodology

3.1 Structural reliability

The study of structural reliability is concerned with calculating and predicting the probability of limit state violation for an engineered system at any stage during its life as defined by Melchers and Beck (2018). Nejad (2018) named the main aim of the structural reliability as an estimation of the failure probability by explicitly taking into account uncertainties of the load, load effect, and resistance. Ultimate limit state reliability in the current work is based on a static safety factor. The important processes used to calculate the ultimate limit state in this paper is illustrated in Fig. 2.

Figure 2 Illustration of the main processes used to calculate the static safety factor, with images extracted from He et al. (2024), Harris et al. (2009), Rezaei et al. (2023), and Solano-Alvarez (2015).

These processes started with simulating the wind turbine in different turbulence wind conditions. From the turbine responses, the loads on the blade bearing are derived. Different dimensions are studied, and material probability is considered. The contact area is calculated, and the safety factor is derived.

Figure 3 Flowchart of reliability assessment in a blade bearing for ultimate limit state.

3.2 Safety factor and failure function

A safety factor is a measure used in engineering and design to provide a margin of safety for structures, materials, or systems under expected loads or conditions. It accounts for uncertainties in the design process, such as variations in material properties, manufacturing tolerances, unexpected loads, and potential degradation over time. In this study, the formula that is used for the failure function has the same form as the static calculation function and static safety factors.

For the reliability model, the static limit state is formulated with the Hertzian contact stress criterion. The approach considers that bearing static damage occurs when the maximum Hertz contact stress exceeds the allowable Hertz contact stress. The DG03 uses a comparison of the maximum contact stress, σ_max, in the limit load condition to the maximum allowable stress of 4200 MPa to define the static safety factor (S₀) (Harris et al., 2009; Stammler et al., 2024).

$$\begin{array}{}\text{(1)}& S_{\mathrm{0}}={\left({\displaystyle \frac{\mathrm{4200}}{{\mathit{\sigma }}_{\text{max}}}}\right)}^{\mathrm{3}}\end{array}$$

Here the maximum contact stress, σ_max, is also expressed in megapascals (MPa), calculated as

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\text{max}}={\displaystyle \frac{\mathrm{1.5}{Q}_{\text{max}}}{\mathit{\pi }ab}},\end{array}$$

and the static safety factor can be rewritten as

$$\begin{array}{}\text{(3)}& S_{\mathrm{0}}={\left({\displaystyle \frac{\frac{\mathrm{4200}\mathit{\pi }ab}{\mathrm{1.5}}}{Q_{\text{max}}}}\right)}^{\mathrm{3}}.\end{array}$$

In Eq. (2), πab is the contact area, which is an ellipse with a semi-major axis a and semi-minor axis b, and Q_max is the maximum ball force. The maximum ball load is calculated from Stammler et al. (2024):

$$\begin{array}{}\text{(4)}& Q_{\text{max}}=\mathrm{0.55}\left({\displaystyle \frac{\mathrm{2}{F}_{\mathrm{r}}}{Z\cos \mathit{\alpha }}}+{\displaystyle \frac{F_{\mathrm{a}}}{Z\sin \mathit{\alpha }}}+{\displaystyle \frac{\mathrm{4.4}M}{D_{\text{pw}}Z\sin \mathit{\alpha }}}\right),\end{array}$$

where F_r, F_a, and M denote the applied radial, axial, and moment loads, respectively.

In Eq. (3), the numerator and denominator are named R and S, respectively. a and b in the R are functions of the applied maximum ball load Q_max; therefore, R is implicitly a function of S. The failure function, g_x, is defined below, where x values are random variables and n is the static safety ratio.

$$\begin{array}{}\text{(5)}& g_x(R,S)={\left({\displaystyle \frac{R}{S}}\right)}^{\mathrm{3}}\le n\end{array}$$

The static safety ratio, n, determines the boundary of the damage. If the failure function value is equal to or smaller than the static safety ratio, the bearing is in a failure state; otherwise, the bearing is in a safe state. In the work of Keller and Guo (2022), it is recommended that the static safety ratio be greater than 1.5; however, Stammler et al. (2024) noted that it seems reasonable to refer to the limit of 1, as in IEC 61400-1 (2019). In order to compute the actual reliability and allow for comparisons, the static factor ratio equal to 1 is considered. Consequently, g_x will become

$$\begin{array}{}\text{(6)}& g_x(R,S)={\left({\displaystyle \frac{R}{S}}\right)}^{\mathrm{3}}\le \mathrm{1}.\end{array}$$

Taking the cube root of Eq. (6) yields the failure function directly as

$$\begin{array}{}\text{(7)}& g_x(R,S)=R-S\le \mathrm{0}.\end{array}$$

The annual probability of failure, P_f, is then obtained from

$$\begin{array}{}\text{(8)}& P_{\mathrm{f}}=P\left(g_x\right(R,S)\le \mathrm{0})=P(R-S\le \mathrm{0}).\end{array}$$

To compute P_f, different methods can be used, including first-order and second-order reliability methods (FORMs and SORMs), as well as the Monte Carlo method (Ditlevsen and Madsen, 2007). In this study, the Monte Carlo simulation method is employed to estimate the probability of failure. Using the Monte Carlo simulation method, a suitably large sample of typical load configurations is simulated from the probabilistic action model. This load configuration sample gives a corresponding sample of load effects at different points of the bearing, and from this sample, the probability distributions of the load effects can be estimated. By using these probability distributions, extreme value studies can subsequently be conducted. In this sense, the probabilistic model uses typical load configurations in its solution procedure rather than difficult choices of “extreme” load configurations as used in the deterministic model (Ditlevsen and Madsen, 2007).

The randomness of the failure function arises from different aspects. Uncertainties in material, forces, and models are some of the main reasons. These sources of randomness can appear in the R representation of the load-capacity term and the S representation of the applied maximum ball load. The randomness in R and S in the current work originated from uncertainty in the materials, dimensions, wind turbulence intensity, and simulation model in ball forces. Figure 3 presents a systematic approach for the static overload reliability assessment of a blade bearing.

The dashed box represents uncertainty in loads that have distinct uncertainty sources – (i) external wind variability and (ii) epistemic uncertainty arising from model, measurement, and stochastic uncertainties – which are sampled separately and then combined to generate the stochastic wind-induced loads. The uncertainties in the measurement and statistics are not considered in this study. Uncertainty in material leads to uncertainty in resistance, while uncertainty in dimension leads to uncertainty in both stress and resistance. The formula for resistance, extracted from Eq. (9) with consideration of uncertainty in material and dimension, is presented below:

$$\begin{array}{}\text{(9)}& R=\left({\displaystyle \frac{\mathrm{4200}{\mathit{\chi }}_{\mathrm{m}}\mathit{\pi }a\left({\mathit{\chi }}_{\mathrm{d}}\right)b\left({\mathit{\chi }}_{\mathrm{d}}\right)}{\mathrm{1.5}}}\right),\end{array}$$

where χ_m and χ_d are uncertainties in material and dimension, respectively. The uncertainties in blade bearing design dimensions, such as pitch diameter, ball diameter, contact angle, and groove conformity, affect the dimensions of the contact area, a and b. It should be noted that a and b are, in addition to the uncertainty in dimension, functions of uncertainty in the loads and the maximum ball forces Q_max. In every Monte Carlo realization, therefore, Q_max is computed first, followed by the evaluation of a and b and then R, ensuring that the dependency R(S) is fully captured. The formula that represents S is presented below:

$$\begin{array}{}\text{(10)}& S=\left(Q_{\text{max}}{\mathit{\chi }}_{\mathrm{f}}\right),\end{array}$$

where χ_f is the uncertainty in the maximum ball force. Q_max contains external loads and bearing dimensions, which are in fact uncertain parameters. Therefore, in its nature, it consists of uncertainty.

The probability of failure including all uncertainties will then be

$$\begin{array}{}\text{(11)}& P_{\mathrm{f}}=P\left(\left({\displaystyle \frac{\mathrm{4200}{\mathit{\chi }}_{\mathrm{m}}\mathit{\pi }a\left({\mathit{\chi }}_{\mathrm{d}}\right)b\left({\mathit{\chi }}_{\mathrm{d}}\right)}{\mathrm{1.5}}}\right)-\left(Q_{\text{max}}{\mathit{\chi }}_{\mathrm{f}}\right)\le \mathrm{0}\right).\end{array}$$

Table 5 PDFs of nominated distribution function (Alavi et al., 2016; Shi et al., 2021; Yu et al., 2019).

3.2.1 Uncertainty in material

The strength of bearings is not a deterministic value. ISO 76 (2006) recommends a Hertz contact stress of 4200 MPa for ball bearings, which is equivalent to a total permanent deformation of 0.0001 of the rolling element diameter. Assessing the uncertainty in the material's strength requires extensive testing. Shimizu et al. (2010) presented the mechanical properties of the AISI 52100 bearing steel. They showed that tensile strength and Rockwell C hardness have a normal distribution with a standard deviation of 5.7 % and 1.1 %, respectively. Lai et al. (2009) presented a model for plastic indentation, and they tested it on 42CrMo4 steel. Their model predicted that the contact pressure causing a plastic indentation of 10⁻⁴D in the through-hardened raceway is 4260 MPa and showed good validation results. In their extended work, Lai (2011) predicted the contact pressure to be 4270 MPa. Lewis et al. (2015) assessed 9 different sources, including 7 bearing manufacturers, and allowable peak Hertzian pressures of 4270 and 3962 MPa for SAE 52100 steel and AISI 440C steel, respectively, based on the mean minimum hardness. Wang and Zhang (2022) performed a reliability analysis on an angular contact ball bearing. They considered allowable yield stress as a strength parameter with a 5 % standard deviation. In another work, the density and modulus of elasticity in the study by Cheng et al. (2020) are considered to follow a normal distribution with a 5 % standard deviation. In their sensitivity analysis, they concluded that the material has the highest reliability sensitivity. Imdad et al. (2024) studied 42CrMo4 steel that was subjected to different heat treatments and hardness levels. The hardness level in different heat treatments has a deviation between 2 % and 6.5 %. The current work considered a normal distribution with a standard deviation of 5.7 % for the uncertainty of material, χ_m.

3.2.2 Uncertainty in dimensions

Blade bearing dimensions affect the loads and resistance in the failure function. Wang and Zhang (2022) performed a sensitivity analysis on four parameters of the bearing, including the ball pitch diameter, ball diameter, and inner and outer raceway groove curvatures with normal distributions and 0.5 % standard deviations. They concluded that, regarding the bearing geometry, the ball diameter has the highest effect on reliability to prevent plastic deformation. In the sensitivity analysis of Cheng et al. (2020) regarding an angular contact ball bearing, the free contact angle and the inner and outer raceway curvatures were treated as random parameters, along with the ball diameter and groove curvatures. All of their random parameters had 0.5 % standard deviations. In a reliability and sensitivity analysis of spherical roller bearings, Wang et al. (2020) considered 4 %, 2.2 %, and 0.5 % standard deviations for roller diameter, pitch circle diameter, and radial clearance, respectively. The effect of dimensions on the reliability of the bearing is presented in the “Sensitivity analysis” section. It should be noted that each dimension was analyzed independently by the relevant distribution.

Table 6 Uncertainty distributions.

3.2.3 Uncertainty in loads

Wang and Zhang (2022), in their work, considered normal distributions with a 5 % standard deviation for axial and radial forces. The same consideration was also seen in the work of Cheng et al. (2020). In contrast, Wang et al. (2020) considered a standard deviation of 2.5 % for the radial load.

The uncertainty in the load at the blade bearing originates mainly from turbulence acting on the wind turbine, and the turbulence has a great contribution to both the safety factor and the fatigue life of the blade bearing, as shown in the work by Rezaei et al. (2023). Different realizations of turbulence, called the “seed number”, produce a Gaussian distribution of TI in the longitudinal wind component due to spatial coherence (Jonkman, 2009).

Each simulation with a specific seed number leads to a time series of distributions of the loads in the balls, while the extreme ball load can be obtained from these series. Different random seed number simulations result in a series of extreme ball loads in the blade bearing, while the turbulence intensity is constant. These extreme loads form a probability distribution function. It is important to assign a proper probability distribution function to these extreme loads.

Figure 4 Annual probability density function of different distribution functions in 15 seeds (left) and 3000 seeds (right).

Figure 5 The exceedance probabilities of different distribution functions in 15 seeds (left) and 3000 seeds (right).

In this study, the following probability distribution functions were considered: the generalized extreme value, gamma, inverse Gaussian, kernel, lognormal, Nakagami, Rician, normal (Nor), the extreme value (EV), and Weibull distributions. In Table 5, the probability density function (PDF) of the nominated distribution function is presented, where x is a random variable. More information on the equations and parameter definitions of generalized extreme value, gamma, kernel, and Weibull can be found in Shi et al. (2021). The parameters of inverse Gaussian, lognormal, and Nakagami are referred to in Alavi et al. (2016). Rician parameters are referred to in Yu et al. (2019).

Table 7 Coefficient of efficiencies of the nominated probability distribution functions in different seed numbers.

In this study, the parameters were calculated by the maximum likelihood estimator (Bain and Antle, 1967) using MATLAB software, and to assess the performance and goodness of fit (GoF) of the distribution functions, the coefficient of efficiency (CE) method was applied. CE is intended to range from 0 to 1, but negative scores are also permitted. The maximum positive score of 1 represents a perfect model; a value of 0 indicates that the model is no better than a one-parameter “no-knowledge” model, in which the forecast is the mean of the observed series at all time steps; negative scores are unbounded; and a negative value indicates that the model performs worse than a “no-knowledge” model (Dawson et al., 2007). The CE indicator is 1 minus the ratio of the sum square error to the statistical variance of the observed dataset about the mean of the observed dataset.

$$\begin{array}{}\text{(12)}& \text{CE}=\mathrm{1}-{\displaystyle \frac{\sum (Q_i-\widehat{Q_i}{)}^{\mathrm{2}}}{\sum (Q_i-\overline{Q}{)}^{\mathrm{2}}}}\end{array}$$

Q_i is the observed data at level i, $\widehat{Q_i}$ is the estimated data at level i, and $\overline{Q}$ is the mean of the observed data.

Different seed numbers were studied to create a distribution function. According to IEC 61400-1 (2019), in ultimate strength analysis, 15 different simulations are necessary for each wind speed from ($V_{\mathrm{r}}-\mathrm{2}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) to cut-out, and six simulations are necessary for each wind speed below ($V_{\mathrm{r}}-\mathrm{2}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). V_r is the rated wind speed that is 11.4 m s⁻¹ according to Table 1. However, to generate coherent turbulent structures, using more than 30 different random seeds for a specific set of boundary conditions is recommended by Jonkman (2009). This study covers a wide range of seed numbers. The generalized extreme value distribution is selected to model the maximum load distribution. The reason behind this is presented in the Results section.

Figure 6 Changes in mean of generalized extreme value distribution due to seed numbers.

3.2.4 Uncertainty in the maximum ball force

The uncertainty in the maximum ball forces arises from the distribution of the forces inside bearings. The flexibility of the bearings and connecting components, hub and blade, can affect the load distributions inside the bearing, as was shown in the works of Menck et al. (2020) and Rezaei et al. (2024). In addition, the results of the maximum ball force equation are not necessarily conservative. It can overestimate or underestimate the actual loads (Stammler et al., 2024).

Menck et al. (2020) developed a finite-element model (FEM) and calculated the bearing load distributions. Moreover, Graßmann et al. (2023) validated the finite-element model with extensive experimental data on the blade bearing. Furthermore, Rezaei et al. (2024) compared the load distribution from the finite-element model and multi-body simulations (MBSs) for the NREL 5 MW and IWES 7.5 MW wind turbines. The average of the maximum ball force differences between the MBS and FEM was 10.8 %. In another work by Leupold et al. (2021), load distributions inside the bearing for two different conditions of finite-element and multi-body simulations were studied, and the average error was 6.5 %.

The maximum ball forces and load distributions in the work by Rezaei et al. (2024) were recalculated with the maximum ball force equation. The differences between the maximum forces from the FEM and the maximum ball equation have a mean and standard deviation of 13 % and 6.5 %, respectively, which is considered an uncertainty in the maximum ball force, χ_f.

The distributions and the mean and standard deviation of model uncertainties based on the above discussions are summarized in Table 6.

Figure 7 P_f with a different (a) ball diameter, (b) pitch circle diameter, (c) raceway conformity, and (d) contact angle. (e) Summary of (a) to (d).

3.3 Description of DLC

It is observed that design load case (DLC) 1.3 of IEC 61400-1, around rated wind speed, has the largest effect on the load of the blade bearings (Rezaei et al., 2023); therefore, it is considered a nominated load case.

DLC 1.3 uses an extreme turbulence model. The DLC covered mean wind speeds from 10 to 13 m s⁻¹ with an interval of 0.5 m s⁻¹. The simulations last 700 s, and the results of the first hundred seconds are not considered. In all DLCs, wind shear is considered according to the sites or related standard conditions.

4 Results and discussion

4.1 Probability distribution function

By increasing the seed number, a bigger number of realizations is created, and the accuracy of the result is higher, but on the other hand, the simulation time will increase. Seed numbers from 15 to 3000 were studied under onshore wind conditions in IEC category IA. The results of the probability density function and exceedance probability of the fitted distributions for border seed numbers 15 and 3000 are depicted in Figs. 4 and 5.

CE indicator results in different seed numbers with nominated probability distribution functions are presented in Table 7.

Figure 8 ULS probability of failure of IEC class I wind conditions in different sample numbers in the Monte Carlo simulation.

The results show that the kernel and generalized extreme value distributions perform best in modeling the extreme load distribution. This dominance started with 150 seeds. The kernel estimator does not have a closed formula, and the generalized extreme value is considered a distribution for modeling the effect of seed number on the distribution of extreme loads. The results also indicate that a low number of seed values cannot accurately represent the true variety of the extreme loads.

Changes in the mean of the GEV function due to the seed number are plotted in Fig. 6.

Figure 9 ULS annual probability of failure of IEC wind conditions.

Figure 10 Load index of IEC wind conditions.

It is very time-consuming to consider the high number of random seeds in the simulations. In addition, changes in the GEV parameters in 300 seeds or more are less than 3 %. A seed number of 300 is the value that is considered for the random seed number in the rest of the results.

4.2 Sensitivity analysis

The probability of failure in the bearing with variation in the ball diameter, pitch circle diameter, conformity, and contact angle is studied. The onshore wind field with an extreme turbulence intensity grade of IA according to IEC 61400-1, DLC 1.3, is considered. A total of 10⁸ samples were considered in the simulation with the Monte Carlo method, and this process was repeated 20 times, forming a cluster.

4.2.1 Ball diameter

The nominal size of the ball diameter is 75 mm. Although the balls are usually manufactured and sorted in a batch with fine tolerances in diameters, it is assumed that the ball diameter can change from fine to very coarse machining, according to ISO 2768-1 (1989). In every analysis, the balls' diameters are assumed to be the same, and a range of diameters was studied. However, the extreme tolerances are not realistic; they can help to observe the trend of changes in reliability. The assumption leads to a 0.15 to 1.5 mm variation in the ball diameter. The failure probability of the bearing with different ball diameters is shown in Fig. 7a. The vertical error bar refers to the maximum and minimum values in each cluster. Increasing the ball diameter increases the reliability of the bearing, which is in good agreement with the results from Wang and Zhang (2022); however, the reliability decreases more sharply in their analysis.

4.2.2 Pitch circle diameter

The nominal size of the pitch circle diameter is 3558 mm. It is assumed that the pitch circle diameter can change from medium (2 mm) to very coarse (8 mm) machining, according to ISO 2768-1 (1989). In order to observe the wider range of diameters, diameters of 3540 and 3576 mm were added to the study. The failure probability of the bearing with different pitch circle diameters is shown in Fig. 7b. The results show that the pitch circle diameter does not have a significant effect on P_f.

Figure 11 ULS annual probability of failure of nominated wind sites.

4.2.3 Raceway conformity

Raceway conformity is the dimensional relationship between the radius of the raceway and the diameter of the ball. The nominal size of the conformity in the current study is 0.53. Bearing manufacturers recommend a value between 0.510 and 0.543 for this ratio (Daidié et al., 2008). In this study, the conformity between 0.515 and 0.545 is studied. The failure probability of the bearing with different raceway conformities is shown in Fig. 7c. The results of P_f show that, with an increase in conformity, reliability decreases sharply. Wang and Zhang (2022) reached the same conclusion, but the decrease in reliability was not as sharp as the results presented. Wang et al. (2016) obtained similar results regarding the maximum Hertzian contact stress in their study on angular contact ball bearings. In order to understand how much the manufacturing of the ball and raceway can affect the reliability of the bearing, it is assumed that the ball and raceways have a fine degree of manufacturing according to ISO 2768-1 (1989), where in our study the tolerance would be 0.15 mm, and the extreme values for conformity would be 0.527 and 0.533. The extreme values for the raceway conformity by this assumption are shown with vertical lines in Fig. 7c. It is observed that raceway conformity has a dominant effect on blade bearing damage in ULS. The result shows that with small changes in groove conformity, the probability of failure increases or decreases significantly. Consequently, the uncertainty in the raceway conformity with normal distribution with a standard deviation of 0.5 % for the uncertainty of the dimension, χ_d, is considered.

Figure 12 Load index of wind sites.

4.2.4 Contact angle

The nominal size of the initial contact angle is 45°. The initial contact angle in this study is referred to as the nominal contact angle, which is in a load-free condition. The contact angle from 25 to 65° is studied. The failure probability of the bearing with different contact angles is shown in Fig. 7d.

The results show that the probability of failure decreases as the contact angle increases. Cheng et al. (2020) obtained the same results for shear stress in their research on angular contact ball bearings.

The initial contact angle, pitch circle diameter, and ball diameter do not have a significant effect on failure probability and are not regarded in the paper as an uncertainty variable with respect to dimension.

4.3 IEC wind conditions

The IEC wind categories I, II, and III in the turbulence intensities of A, B, and C under onshore and offshore conditions are studied. Simulations with different numbers of samples were performed. The result of P_f in IEC class I is depicted in Fig. 8.

As shown by the results, P_f converges after 10⁷ samples in all wind conditions. To account for a wide range of samples, 10⁸ samples were considered in the simulation with the Monte Carlo method, which was repeated 20 times (20 clusters of 10⁸ samples). P_f is the average of 20 clusters. The IEC wind condition probability of failure based on the reference parameters of the bearing stated in Table 2 is illustrated in Fig. 9.

IEC onshore A class IA has the highest probability of failure, and the overall probabilities of failures of IEC wind configurations are of the order of 10⁻⁵. By increasing the annual mean wind speed, reliability decreases; however, turbulence intensity has a more significant effect, and reliability decreases when the turbulence intensity increases. While the P_f mean value varies between $\mathrm{2.09}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{4.95}\times {\mathrm{10}}^{-\mathrm{6}}$, the standard deviation of the clusters varies between $\mathrm{5.56}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{7}}$. The variance of the results is too small, and it indicates that the clusters are closer together, suggesting less diversity and more consistency. IEC 61400-1 (2019) and IEC 61400-8 (2024) set the target value for the nominal failure probability for structural design, for both extreme and fatigue failure modes, over a reference period of 1 year, at $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ for component class 2. Component class 2 is a “safe-life” structural component whose failure may lead to the failure of a major part of a wind turbine (IEC 61400-8, 2024). All the wind configurations have a lower failure probability than the target value.

IEC 61400-1 (2019) recommends 15 seed numbers in the ultimate analysis. It is shown that 15 seed numbers cannot represent the behavior of the probability distribution of the extreme loads. To investigate further, the load index is introduced. The load index is the ratio of the extreme ball loads in 300 seed numbers to the extreme ball loads in 15 seed numbers. The load index results for the IEC wind categories are plotted in Fig. 10. The results show that the extreme load calculation with 15 seeds has an error between 2 % and 11 %. This exercise is referred to as a code–site comparison.

4.4 Wind sites

The wind sites that were introduced previously were studied. The annual probability of failure for the nominated wind sites is illustrated in Fig. 11. For comparison, the maximum P_f value among the IEC wind conditions is also included in the figure.

The results show that most wind sites exhibit a higher probability of failure in the blade bearing under ULS conditions than indicated by the IEC categories. The reliability at the Sujawal and Aysha wind sites is far lower than that of the IEC standards. In addition, these two sites have higher failure probabilities than the failure target value for component class 2 (IEC 61400-8, 2024). The standard deviation of the clusters varies between $\mathrm{4.82}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{2.33}\times {\mathrm{10}}^{-\mathrm{7}}$, and the results of the clusters are consistent. These two wind sites have annual wind speeds between 7.5 and 8.5 m s⁻¹ and are categorized into the IEC II class, while their P_f is higher than that of the IEC I class wind sites. These high P_f values are the result of high turbulence intensity, as addressed in Table 4. The load index results for the wind sites are plotted in Fig. 12.

The results show that the extreme load calculation with 15 seeds has an error between 3 % and 30 %. This observation has important implications for reliability analysis and simulation practice. It shows that using only 15 seeds, as typically done in standard simulations, can lead to a non-negligible underestimation of extreme ball loads, especially in complex wind site conditions. The load index provides a simple and effective measure to quantify this underestimation. It can help practitioners evaluate whether their simulation setup sufficiently captures the load extremes, and when limited seed numbers are used, the load index can offer a basis for applying correction factors to improve the accuracy of failure probability estimates. Additionally, the variation in load index across different wind sites emphasizes the need for site-specific assessment when evaluating blade bearing reliability.

5 Conclusions

The presented work evaluates the probability of static overload in a double-row, four-point contact ball blade bearing under ultimate limit state conditions. The analysis is based on a structural reliability approach, quantifying the likelihood of permanent deformation and the static safety factor, which is the ratio between the maximum permissible Hertzian contact stress and the maximum contact stress. The NREL 5 MW reference wind turbine is considered with the extreme turbulence wind model in the design load case. The Monte Carlo method is used to calculate the probability of failure. It is shown that the generalized extreme value is a suitable distribution to simulate the probability of the extreme ball load distribution. Increasing the number of seeds in the turbulence simulation improves the accuracy of the estimated probability distribution. It is observed that by considering 15 seed numbers, as proposed in the standards and guidelines, the distribution of the loads is not represented.

The probability of blade bearing failure regarding variations in four main dimensions was studied. In this regard, raceway conformity has the highest contribution to the probability of static overload of the blade bearing, followed by ball diameter. Strictly controlling these parameters can lead to higher reliability in the bearing regarding the ultimate limit state.

The P_f results of different IEC wind conditions show that the IEC IA onshore category has the highest probability of failure, and the IEC IIIA offshore category has the lowest. However, the wind class (I, II, III), whether onshore or offshore, has minimal impact on the probability of failure, while turbulence intensity (A, B, C) has a significant effect on reliability. The probabilities of failure for the selected onshore and offshore wind sites are generally higher than those of IEC sites. This indicates that IEC-designed turbines may experience increased risk of static overload in blade bearings when applied to certain real wind sites. It reinforces the need for site-specific assessment considering turbulence and manufacturing variability. It also shows the necessity of assessing the blade bearing in every wind site condition according to extreme turbulence wind.