Semi-analytical methodology for fretting wear evaluation of unlubricated pitch bearing raceways under operative and non-operative periods

Current methodologies for the evaluation of fretting in pitch bearing raceways only consider damage induced by oscillating control movements. However, pitch bearings can remain static for long operative and non-operative periods, where load fluctuations cause fretting damage due to the small movements and deformations at the contact. In this work a semi-analytical methodology is presented and an analysis of the fretting phenomenon in pitch bearing raceways in the absence of lubricant under both productive and non-productive periods of static pitch control is accomplished. The analysis is performed following the energy-based wear model and the calculation of a total of 30 times series of 10 min of duration concerning winds speeds from 3 to 25 m s⁻¹. As a result, critical locations based on the dissipated energy are identified for the different wind speeds and the contributions of operational and non-operational times, as well as the prediction of damage shapes on the raceway for both cases, are evaluated independently and compared.

1 Introduction

The wind industry has grown remarkably over recent decades compared to the general industry due to the high demand for clean and renewable energy (Nematollahi et al., 2016; Rodríguez et al., 2019). This fact has been accompanied by new advances in design, materials, and manufacturing techniques that have ultimately resulted in wind turbines of larger dimensions and, consequently, larger components (Serrano et al., 2016). This fact and extreme operational conditions have led to complex scenarios where the methodologies for common components are not reliable. While new calculation procedures and design guidelines have been proposed for the analysis of the reliability of the pitch bearing raceway in terms of rolling contact fatigue and fretting damage, these techniques are supported by simple methods or high safety factors that generally result in conservative calculus and low design optimization (Harris et al., 2009; Stammler et al., 2024; Portugal et al., 2017; Houpert, 1999; Schwack et al., 2016; Heras et al., 2017; Olave et al., 2010; Lopez et al., 2019). This fact has promoted the usage of experiments to test the reliability of the wind turbine components under their different failure modes (Stammler and Geibel, 2024; Stammler et al., 2019; Menck et al., 2020; Schwack et al., 2021a). This technique has been applied to real-scale components, and therefore, it has provided effective and reliable results. However, it involves high costs due to the need for real-scale components, materials, and equipment, and these costs are even increasing due to continuous component growth. This raises the need for design methodologies and scaled testing approaches whose results are valid for real-scale applications (Olave et al., 2019).

Pitch bearings are typically large four-point contact ball (4PCB) bearings; they link the blades to the hub and allow pitch control to locate the blade in the optimal position to maximize lift of the airfoils and to reduce the structural loads when the wind is too high (Bossanyi, 2003). Due to the pitch control strategy, these bearings usually oscillate, and, eventually, the amplitude of the movement can be low enough (Harris et al., 2009) to result in rotational fretting at the bearing raceway (Schwack et al., 2016; Wandel et al., 2023; Stammler et al., 2019) as a consequence of the reciprocating rolling motion of the balls (Cai et al., 2020) and the evacuation of the lubricant out of the contact between the rolling element and the raceway. The fretting damage is also known as standstill marks (small amplitudes, $x/\mathrm{2}b<\mathrm{1}$) or false brinelling (small amplitudes, $x/\mathrm{2}b>\mathrm{1}$) (Presilla et al., 2023). This phenomenon has been extensively studied for different authors, whereby analytical (Cubillas et al., 2021b), numerical (Schwack et al., 2018; Lin et al., 2022) and experimental (Grebe et al., 2011, 2020; Schwack et al., 2021b) (Pittroff, 1961; Stammler, 2023) methods have been used to study the effects of these control movements on the fretting damage. As a result, it has been concluded that oscillatory movements of enough amplitude promote the relubrication of the raceway and aid in avoiding fretting damage. Therefore, the premeditated use of the control movements, called protection runs, has been experimentally studied to observe the movements' effects on the damage, with the aim of avoiding damage (Stammler et al., 2019).

Under some circumstances, such as pre-commissioning, low wind speeds, or safe operating stops, pitch control can remain static for long periods of time. Under such a scenario, pitch bearing must still accommodate tilting moments exerted by the wind and the blades' weight that are transmitted between the rings through the balls. As a consequence, balls experience variable loads, which may cause radial fretting (Zhu et al., 2006; Cubillas et al., 2021a), and additionally, this variable load causes small deformations at the ball and the rings' contact that result in small rolling movements (Olave et al., 2010) in the transversal direction of the raceway, squeezing out the lubricant, and, ultimately, they may cause rotational fretting.

Despite the efforts invested in the literature into the analysis of the pitch control moments on pitch bearing raceways, the effects of the variable loads when the pitch control remains static have not been studied. In our previous work, we developed and validated formulations for the analysis of radial fretting (Cubillas et al., 2021a), rotational fretting (Cubillas et al., 2021b), and the combination of both (Cubillas et al., 2021c) in small bearings.

The objective of this work is to propose a methodology to analyse the relative contributions of productive and non-productive periods of wind turbine operation to pitch bearing fretting damage, including identification of the most likely location of damage on each raceway. For this task, the 5 MW NREL reference wind turbine has been taken as a case study, and a total of 30 time series of different design load cases (DLCs) of normal production (DLC 1.2) and non-productive condition periods (DLC 6.4) from IEC 61400-3-1 have been evaluated through an energy-based wear model.

In this work, the damage evaluation has been addressed using energy-based wear models that have previously demonstrated adequate correlation to damage (Brinji et al., 2020; Schwack et al., 2018; Cubillas et al., 2021b). Although the present analysis does not consider the effect of the lubricant, pitch bearings are likely to be in boundary lubrication conditions when the pitch angle is not changing in non-operative periods and even in some low-wind-speed operations in which the model is valid. Therefore, the authors believe that the proposed method provides a reliable framework for the prediction of the most critical areas on the raceway where fretting damage can occur, as well as providing a basis for comparison between the severity of different wind speeds and operating conditions.

2 Methodology

In previous work (Cubillas et al., 2021c), the authors developed and validated a complete methodology to predict fretting damage in static angular bearings subject to variable loads. In this work, this methodology is extended and adapted to the 4PCB bearing problem.

Figure 1 shows the flowchart of the complete methodology. The required inputs are the wind time series containing the value of the bearing reaction along the time M_x(t), M_y(t), and the geometry. The maximum value of the bearing reaction and the geometry are used to perform a series of experiments with finite-element model (FEM) simulations that determine a surface response with the values of every ball reaction and contact angle as function of the bearing reaction, ball ID ($i=\mathrm{1}:z$), row ($j=\mathrm{1}:\mathrm{2}$), and contact ID ($k=\mathrm{1}:\mathrm{4}$); see Fig. 2.

Figure 1 Flowchart of the proposed semi-analytical formulation showing the necessary calculation process with the loops involved to calculate the fretting damage, where the suffix i defines the number ID of the ball, the suffix j defines the row ID, the suffix k defines the contact location, and the suffix t defines the time step. Q is the normal load over a contact, α is the contact angle, r_c is the contact-deformed radius, r is the distance from a contact point to the rolling axis, r^′ is the effective rolling radius, W is the local rigid slip, T is the tangential stresses, S is the local sliding, and E_r is the rolling energy.

Figure 2 Generic 4PCB bearing of pitch radius r_pw; ball radius r_w; and the definition of the contact through the suffix j, which defines the ID of the row, and the suffix k, which defines the ID of the contact location.

Once the response surface is created, it is possible to calculate the contact reaction, Q(t), and the contact angle, α(t), of every contact point of a ball in a specific time, and then the contact reaction provides enough information to calculate the contact deformations; the variation in the contact angle allows us to determine the ball motion. The kinematics of the ball are calculated following the authors' previous work (Cubillas et al., 2021c), with an iterative solver for rolling-friction minimization. This is implemented using the Nelder–Mead simplex algorithm, as described in Lagarias et al. (1998), with a tolerance applied to the function value of 10⁻⁴. Within this iterative solver, the local rigid slip, tangential stresses, and sliding are calculated. Finally, the fretting damage and the probability of damage initiation are calculated. Each of these steps is described in the following subsections.

2.1 Wind turbine: 5 MW NREL reference turbine

For the purpose of analysis, the 5 MW NREL case study has been selected (Jonkman et al., 2009). This model is a utility-scale multi-megawatt turbine, conventionally three-bladed and variable blade-pitch-to-feather-controlled. It has been used as a reference turbine by research teams throughout the world to standardize baseline offshore wind turbine specifications and to quantify the benefits of advanced land- and sea-based wind energy technologies (Zuheir et al., 2019; Cherubini et al., 2021; Halawa et al., 2018).

The bearing case study is a double-row 4PCB bearing attached to the hub and the blade through bolts. Table 1 summarizes the values of the main parameters of the bearing.

Table 1 Bearing geometrical characteristics.

2.2 Pitch control and time series

Pitch control is usually collective pitch control (CPC) or individual pitch control (IPC) (Lopez et al., 2019). While CPC defines the position of the blades simultaneously, the IPC defines the position of each blade individually. Both controls result in oscillating movements, and besides this classification, the final routines of oscillations can be as specific as the turbine design, the world location, and the climate conditions. For this work, a reference control with five seeds of wind time series of the normal turbulence model (NTM) from winds of 3 to 25 m s⁻¹ according to the design load cases (61400-3, 2009) DLC 1.2 (normal production) and DLC 6.4 (parked turbine with idling rotor) is available with an individual duration of 10 min at 20 Hz. Figure 3 shows the evolution of the power production with the wind speed and the value of the pitch bearing angle rate where the pitch control remains static for values under 7 m s⁻¹.

Figure 3 Pitch control performance of the 5 MW NREL reference wind turbine as a function of the wind to a maximum 25 m s⁻¹ speed and the measured power (Lopez et al., 2019).

Following the above observations, only series under 7 m s⁻¹ are considered for DLC 1.2, series of 3, 5, and 7 m s⁻¹, and all available series are considered for DLC 6.4, series of 3, 11, and 25 m s⁻¹. The aeroelastic wind time series are decomposed into bearing reactions using the software Bladed (DNVGL-ST-0437, 2016).

2.3 Ball load distribution and contact angle problem

As described in Cubillas et al. (2021c), an FEM is necessary to calculate the contact reaction, Q, and the contact angle, α, of every contact of every ball under a specific bearing load scenario.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }Q=f\left(M_x,M_y,i,j,k\right)\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=f\left(M_x,M_y,i,j,k\right)\end{array}$$

However, this method is time-consuming, and the calculation of every time increment of the time series is not practical. To manage this inconvenience, an experiment design is accomplished covering all possible load cases (see Table 2) that allows for the creation of a response surface. Radial forces and axial forces are neglected in this analysis as they have been demonstrated not to have a considerable effect on the bearing load distribution (Portugal et al., 2017). A total of 8 simulations are performed where the load is gradually applied in 10 load increments, resulting in a total of 80 load cases.

Table 2 Design of experiments (DoE) for the development of a surface response of the contact reaction and contact angle.

The FEM is built in Ansys considering not only the bearing geometry but also the geometry of the blade, bolts, and stiffeners. Nevertheless, symmetric behaviour is assumed, and only one blade and one-third of the hub is considered; see Fig. 4a. It is important to mention that the results in the areas around the circumference may be affected by this assumption by not considering the cross-effects at the root of the blade. The bearing and stiffeners are structural steel; the hub is cast iron; and, in the absence of a realistic blade design to carry out the calculation, a laminate circular tube with a constant section equal to the bearing ring has been incorporated to apply a geometric offset where the force can be applied. Table 3 summarizes the elastic properties of the materials mentioned.

Figure 4 FEM of the 5 MW NREL reference wind turbine formed by the hub, the pitch bearing, the blade, and the inner and outer stiffeners: (a) simplification of the original geometry considering symmetric behaviour of the hub; (b) application of the tilting moments through a remote point located at the bearing centre; and (c) detail of the bearing mesh and the simplifications of balls and bolts by a spring with non-linear and linear behaviour, respectively.

Table 3 Elastic properties of the wind turbine components.

The entire mesh is composed of second-order elements, and a total of 1 562 761 elements are used in the entire model, with a maximum value of the element size of 10 mm in the bearing. Balls are replaced by two linear elements with variable stiffness connected to both bearing raceways (Daidié et al., 2008). To achieve this aim, rectangular partitions along the bearing raceway of dimensions (a_max, b_max) are used. In these partitions, an ordered mesh of 2×8 elements is generated, and this is fixed to the end of the linear element at its opposite point, as seen in Fig. 4c. Analogously, bolts are replaced by linear elements considering the minimum section; see Fig. 4b.

Moment is applied through a remote point located at the bearing centre and linked to the blade with flexible behaviour. In this way, the decomposition of the moment is prevented, and a true moment is ensured. According to the previous assumption of symmetric behaviour of the hub, frictionless supports are applied to the cut surfaces to prevent normal deformations and allow deformations along the symmetry plane. Frictional contacts are used at the contacts between the bearing rings, stiffeners, and blade, with a 0.3 coefficient of friction with asymmetric behaviour and augmented Lagrange formulation. Bolt extremes are connected to the bearing rings, blade, and hubs with a multi-point constraint contact (MPC), and 608 kN of pre-tension is applied to the bolts based on a 10.9 grade.

2.4 Contact problem

The Brewe and Hamrock formulation (Brewe and Hamrock, 1977) can be used to define the contact area through the major and minor ellipse axis, a and b, respectively (see Fig. 5a):

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }a=a^{*}\cdot {\left({\displaystyle \frac{\mathrm{3}\cdot Q_{ijk}\cdot \left(\mathrm{1}-v^{\mathrm{2}}\right)}{S_{ijk}\cdot E}}\right)}^{\mathrm{1}/\mathrm{3}},\text{(4)}& {\displaystyle }& {\displaystyle }b=b^{*}\cdot {\left({\displaystyle \frac{\mathrm{3}\cdot Q_{ijk}\cdot \left(\mathrm{1}-v^{\mathrm{2}}\right)}{S_{ijk}\cdot E}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where E is the elastic modulus of the material; v is Poisson's ratio; S_ijk is the sum of the contact curvatures of the groove curvature at row j and raceway k; Q_ijk is the normal load at the ball i, row j, and contact k; and a^* and b^* are the dimensionless semiaxis quantity of the major and minor semiaxes, respectively. Appendix A describes the calculation of the above variables.

Figure 5 Contact area and rolling motion of the ball. (a) Definition of the elliptical contact area through the major ellipse semiaxis a, and the minor semiaxis, b, and the local coordinates y^′ aligned with the longitudinal direction of the raceway and x^′ aligned with the transversal direction of the raceway; (b) depiction of the rolling motion of the ball in the raceway transversal direction.

Additionally, both bearing raceway and ball deform to finally coincide in a transversal curvature (Hertz, 1882), r_c:

$$\begin{array}{}\text{(5)}& r_{\mathrm{c}}={\displaystyle \frac{\mathrm{2}\cdot r_{\mathrm{w}}}{\mathrm{1}-s}}.\end{array}$$

2.5 Transversal rolling problem

Following our previous work (Cubillas et al., 2021c), the distance of the ball travelled through the raceway, Δx (see Fig. 5a), can be calculated as

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }x=\mathrm{\Delta }\mathit{\alpha }\cdot r_{\mathrm{c}},\end{array}$$

and the creepage, ε, is calculated as

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_x={\displaystyle \frac{r-r^{\prime }}{r^{\prime }}},\end{array}$$

where r^′ is the effective radius (calculated in the iterative process of energy minimization) and r is the distance from any point to the rolling axis:

$$\begin{array}{}\text{(8)}& r\left(xy\right)=\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-x{{}^{\prime }}^{\mathrm{2}}}-\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}+\sqrt{r_{\mathrm{w}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}.\end{array}$$

Here x^′ and y^′ are the local coordinates of the ball according to Fig. 5a, a(y) is the distance from a point to the ellipse bound in the x direction, and b(x) is the distance from a point to the ellipse bound in the y^′ direction. They can be calculated as

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }a\left(y^{\prime }\right)=\sqrt{a^{\mathrm{2}}\cdot \left(\mathrm{1}-{\displaystyle \frac{y{{}^{\prime }}^{\mathrm{2}}}{b^{\mathrm{2}}}}\right)}\text{and}\text{(10)}& {\displaystyle }& {\displaystyle }b\left(x^{\prime }\right)=\sqrt{b^{\mathrm{2}}\cdot \left(\mathrm{1}-{\displaystyle \frac{x{{}^{\prime }}^{\mathrm{2}}}{a^{\mathrm{2}}}}\right)},\end{array}$$

respectively.

2.6 Tangential friction problem

Following Kalker (1981), the relative motion of two points, S (local slip), can be decomposed in the rigid displacements, W (local rigid slip), and the elastic displacements, U (local elastic slip):

$$\begin{array}{}\text{(11)}& \mathit{S}=\mathit{W}-\mathit{U}.\end{array}$$

Considering an initial condition of no slip, the local elastic slip is equal to the local rigid slip:

$$\begin{array}{}\text{(12)}& \mathit{S}=\mathrm{0}\to \mathit{U}=\mathit{W}.\end{array}$$

The following step is to determine the local rigid motion at the contact points:

$$\begin{array}{}\text{(13)}& \mathit{W}\left(W_x^{\prime },W_y^{\prime }\right)=W_{i,j,k}^{\mathrm{NL}}\left(W_x^{\mathrm{NL}},W_y^{\mathrm{NL}}\right)+W_{i,j,k}^{\mathrm{Ro}}\left(W_x^{\mathrm{Ro}},W_y^{\mathrm{Ro}}\right),\end{array}$$

where W^NL is the local rigid slip caused by the variable normal load and W^Ro is the local rigid slip caused by the rolling effects.

According to our previous work (Cubillas et al., 2021c), the local rigid slip caused by the normal load is calculated as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{W}_x^{\mathrm{NL}}& =r_{\mathrm{c}}\left({\tan }^{-\mathrm{1}}\left({\displaystyle \frac{x^{\prime }}{\sqrt{r_{\mathrm{p}}^{\mathrm{2}}-x{{}^{\prime }}^{\mathrm{2}}}}}\right){\displaystyle \frac{{\tan }^{-\mathrm{1}}\left(\frac{a\left(y^{\prime }\right)}{\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}}\right)}{{\tan }^{-\mathrm{1}}\left(\frac{a\left(y^{\prime }\right)}{\sqrt{r_{\mathrm{p}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}}\right)}}\right.\\ & \left.-{\tan }^{-\mathrm{1}}\left({\displaystyle \frac{x^{\prime }}{\sqrt{r_{\mathrm{b}}^{\mathrm{2}}-x^{\mathrm{2}}}}}\right){\displaystyle \frac{{\tan }^{-\mathrm{1}}\left(\frac{a\left(y^{\prime }\right)}{\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}}\right)}{{\tan }^{-\mathrm{1}}\left(\frac{a\left(y\right)}{\sqrt{r_{\mathrm{w}}^{\mathrm{2}}-a(y^{\prime }{)}^{\mathrm{2}}}}\right)}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{W}_y^{\mathrm{NL}}& =r_{\mathrm{c}}\left({\tan }^{-\mathrm{1}}\left({\displaystyle \frac{y^{\prime }}{\sqrt{r_{\mathrm{b}}^{\mathrm{2}}-y{{}^{\prime }}^{\mathrm{2}}}}}\right){\displaystyle \frac{{\tan }^{-\mathrm{1}}\left(\frac{b\left(x^{\prime }\right)}{\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-b(x^{\prime }{)}^{\mathrm{2}}}}\right)}{{\tan }^{-\mathrm{1}}\left(\frac{b\left(x\right)}{\sqrt{r_{\mathrm{w}}^{\mathrm{2}}-b(x^{\prime }{)}^{\mathrm{2}}}}\right)}}\right.\\ & \left.-y{\displaystyle \frac{{\tan }^{-\mathrm{1}}\left(\frac{b\left(x^{\prime }\right)}{\sqrt{r_{\mathrm{c}}^{\mathrm{2}}-b(x^{\prime }{)}^{\mathrm{2}}}}\right)}{b\left(x\right)}}\right).\end{array}\end{array}$$

On the other hand, the local rigid slip caused by the rolling motion considering the transient phenomena (Al-Bender and De Moerlooze, 2008) can be calculated as follows.

If $a\left(y^{\prime }\right)-\mathrm{\Delta }x\le x$ and $x<a\left(y^{\prime }\right)$,

$$\begin{array}{}\text{(16)}& W{x}_{\mathrm{R}x}=\mathit{\epsilon }\cdot (x^{\prime }-a(y^{\prime }\left)\right).\end{array}$$

If $a\left(y^{\prime }\right)-\mathrm{\Delta }{x}^{\prime }>x^{\prime }$ and $x^{\prime }\ge -a\left(y^{\prime }\right)$,

$$\begin{array}{}\text{(17)}& W{x}_{\mathrm{R}x}=\mathit{\epsilon }\cdot \mathrm{\Delta }{x}^{\prime }.\end{array}$$

2.7 Tangential stresses, sliding, and rolling energy

Kalker's simplified theory (Kalker, 1982) is applied to determine the tangential stress T:

$$\begin{array}{}\text{(18)}& \mathit{T}(x^{\prime },y^{\prime })={\displaystyle \frac{\mathit{W}(x^{\prime },y^{\prime })}{L}},\end{array}$$

where L is the flexibility parameter calculated as

$$\begin{array}{}\text{(19)}& L={\displaystyle \frac{\mathrm{8}\cdot a}{\mathrm{3}\cdot C\cdot G}}.\end{array}$$

Here G is the shear modulus and C is the creepage coefficient. The creepage coefficient is tabulated as a function of Poisson's ratio and the contact area; however, for the numerical implementation of the method, a polynomial regression was taken from the tabulated values:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}C& =-\mathrm{1.60}\cdot {\mathrm{10}}^{-\mathrm{3}}\cdot {\left({\displaystyle \frac{a}{b}}\right)}^{\mathrm{3}}+\mathrm{1.08}\cdot {\mathrm{10}}^{-\mathrm{2}}\cdot {\left({\displaystyle \frac{a}{b}}\right)}^{\mathrm{2}}\\ & +\mathrm{0.89}\cdot \left({\displaystyle \frac{a}{b}}\right)+\mathrm{3.44}.\end{array}\end{array}$$

However, tangential stresses must be reconsidered at those points where slip occurs if $T\left(xy\right)>\mathit{\mu }\cdot P(x,y)$.

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }T(x^{\prime },y^{\prime })=\mathit{\mu }\cdot P(x^{\prime },y^{\prime }),\text{(22)}& {\displaystyle }& {\displaystyle }U(x^{\prime },y^{\prime })=\mathit{\mu }\cdot P(x^{\prime },y^{\prime })\cdot L,\end{array}$$

where μ is the value of the coefficient of friction (CoF) and P is the contact pressure:

$$\begin{array}{}\text{(23)}& P(x^{\prime },y^{\prime })={\displaystyle \frac{\mathrm{3}\cdot Q}{\mathrm{2}\mathit{\pi }\cdot a\cdot b}}\cdot \sqrt{\mathrm{1}-{\left({\displaystyle \frac{x^{\prime }}{a}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{y^{\prime }}{b}}\right)}^{\mathrm{2}}}.\end{array}$$

Then, reconsidering Eq. (11), the local slip is as follows:

$$\begin{array}{}\text{(24)}& S(x,y)=W(x,y)-U(x,y).\end{array}$$

Finally, the rolling energy, E_R, is calculated from the local terms:

$$\begin{array}{}\text{(25)}& E_{\mathrm{R}}=\int \int S\cdot T_{\mathrm{l}}+U\cdot T\mathrm{d}x\mathrm{d}y.\end{array}$$

2.8 Fretting damage indicators

An energy-based model has been proposed and validated (Fouvry et al., 2003) for fretting applications, considering the dissipated friction energy where the wear volume, V_w, is proportional to the total accumulated friction energy (TFE), E_F, and the wear coefficient, k_w:

$$\begin{array}{}\text{(26)}& V_{\mathrm{w}}=k_{\mathrm{w}}\cdot E_{\mathrm{F}},\end{array}$$

where E_F is calculated from local terms as

$$\begin{array}{}\text{(27)}& E_{\mathrm{F}}=\int \int S(x^{\prime },y^{\prime })\cdot T(x^{\prime },y^{\prime })\cdot \mathrm{d}x\mathrm{d}y.\end{array}$$

The dissipated energy can be also expressed as a distribution through the density of the accumulated friction energy (DFE), ρ_F:

$$\begin{array}{}\text{(28)}& {\mathit{\rho }}_{\mathrm{F}}(x^{\prime },y^{\prime })=S(x^{\prime },y^{\prime })\cdot T(x^{\prime },y^{\prime }).\end{array}$$

Then both the dissipated energy, E_F, and its density, ρ_F, are instantaneous variables in local coordinate systems. To achieve the cumulative evolution of damage in the global coordinate system, these are translated into cumulative variables: the total accumulated friction energy, TFE, and maximum density of friction energy, MDFE. Thus, considering the change in position due to the variation in the contact angle,

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle }\mathrm{TFE}=\underset{\mathrm{0}}{\overset{t}{\int }}{E}_{\mathrm{F}}\left(t\right)\cdot \mathrm{d}t=\underset{\mathrm{0}}{\overset{t}{\int }}{\mathit{\rho }}_{\mathrm{F}}\left(t\right)\cdot \mathrm{d}t,\text{(30)}& {\displaystyle }& {\displaystyle }\mathrm{MDFE}=\max \left(\underset{\mathrm{0}}{\overset{t}{\int }}{\mathit{\rho }}_{\mathrm{F}}\left(t\right)\cdot \mathrm{d}t\right).\end{array}$$

In accordance with this work, multiple authors have found analogues result in multiple fretting wear applications and in particular for false brinelling prediction in roller bearings (Fallahnezhad et al., 2018, 2019; Brinji et al., 2021), angular bearings (Schwack et al., 2018), and thrust bearings (Cubillas et al., 2021c, 2022). However, the value of k_w has shown dissimilar results as it is affected by multiple conditions such as frequency, sliding distance, or pressure (Brinji et al., 2021), but it will also be greatly affected by lubrication conditions. Considering this fact and keeping in mind the objective of identifying critical areas and adverse behaviour, the damage assessment is performed through energetic parameters only, TFE and MDFE.

3 Results and discussion

In this section the results from the simulations are presented. First, the load and contact angle distributions are studied to determine the areas that are exposed the most to variable loads and variable contact angles. Then the distribution of the damage is analysed and compared to the load distribution to determine its effects on the final damage and the effects of the wind. Finally, the results of the probability of damage initiation are presented, and a cumulative model for variable winds analysis is proposed.

3.1 Load distribution and contact angles

Figure 6 shows the distribution of the mean and the amplitude load; the mean and the amplitude contact angle; and the damage indicators – the total accumulated friction energy, TFE, and the maximum density of friction energy, MDFE – after the simulation of the 15 wind series of 10 min of DLC 1.2 and DLC 6.4, respectively. For a clearer presentation of the results, contacts are classified according to the load diagonal, the row, and the ring (inner and outer); according to the high number of contact points and the similarity between adjacent locations, the results correspond to the mean value of two adjoining balls with the objective of facilitating the presentation of the data. Furthermore, considering the large difference between the damage values obtained in production and non-productive times, it has been decided to use different scales.

Figure 6 Summary of the results of the analysis of the mean load, amplitude load, mean angle, amplitude angle, the total accumulated friction energy (TFE), and the maximum value of the maximum density of friction energy (MDFE), where results are classified according to the row location and the load diagonal and for the different design load cases, DLC 1.2 and DLC 6.4, moment might amplify the effects of the heterogeneous stiffness, resulting in different locations of the maximum values of the loads and the contact angles.

In this section, we first focus on the analysis of the loads and the contact angles. The results show the presence of variable loads and rolling motion (as consequence of the variable contact angle), and therefore the potential presence of radial and rotational fretting damage. Previous classification into load diagonals, rows, and rings allows us to easily observe the same value of the mean and amplitude load and mean and amplitude contact angle at the inner and the outer rings as the contact reactions are equal in magnitude with opposite directions.

Means and amplitude loads and contact angles show the expected symmetric behaviour with respect to a crossing axis from 150 to 330° (see dashed black lines in Fig. 6). Additionally, the values of the mean load and the mean angle on one side and the amplitude load and the amplitude angle on the other side, show some correlation. However, notable differences can be observed as the maximum values are not always located at the same angular position, and another interesting point is the dissimilar distribution between rows that can be observed through the maximum values at row 1 and row 2 being significantly different. These facts show the important effect of the heterogeneous stiffness of the structure and endorse the usage of the FEM to calculate the bearing reactions, as shown in Olave et al. (2010). The distribution of the mean load with the highest values located at D1 at 150° and D2 at 315° shows the effect of a permanent tilting moment of similar components M_x and M_y attributed to the mean load exerted by wind. The highest values of the mean load are 30 kN for DLC 1.2 and 12 kN for DLC 6.4 at D1R2 and D2R1. Similarly, higher values of the mean angle can be found in nearby locations at 120° in D1 and at 320° in D2, with similar maximum mean values of 58° at D2R1 and D1R2 for DLC 1.2 and 50° for DLC 6.4. On the other hand, the results of the analysis of DLC 6.4 show a more homogeneously distributed load over the bearing, where the lower value of the tilting is.

On the other hand, the higher values of the amplitude load are located at 60, 150, and 225° in D1 and 30 and 270° in D2, with a maximum value of 13 kN for DLC 1.2 and 4 kN for DLC 6.4. This distribution reveals the presence of a complex scenario with an alternative tilting moment of components M_x and M_y caused by the alternative load of the blade weight due to the turbine movement and the fluctuation in the wind. The regions of high variable load reveal critical regions for the occurrence of radial fretting. Similarly, the higher values of the amplitude of the contact angle can be found at 90 and 225° at D1 and 20 and 300° at D2, with a maximum value of 6° for DLC 1.2 and 1.9° for DLC 6.4 at D2R1. The amplitude of the contact angle indicates the presence of rolling motion of the balls, and consequently, these regions are critical for the development of rotational fretting.

3.2 Damage indicators

In Fig. 6, both damage indicators, the total accumulated friction energy, TFE, and maximum density of friction energy, MDFE, show similar distributions of the damage, with the TFE maximum at 60, 135, and 180° at D1 and at 30 and 275° at D2 and MDFE at 275° in D2R1 for DLC 1.2 and at 290° in D2R1 for DLC 6.4. These damage distributions retain certain similarities with the mean and amplitude load distributions and the mean and amplitude angle pointed out in the previous section. However, the correlation of the distribution of damage seems to have the best fit with the distribution of the amplitude load. Therefore, results might indicate that the main source of damage is the radial fretting caused by the variable load of the wind and the blade weight dynamics. This fact is aligned with previous work (Cubillas et al., 2021b), where the analysis of an angular bearing indicated that 70 % of the fretting damage caused for a variable load was produced by radial fretting. In addition, TFE and the value of MDFE show damage distributions where the inner ring is significantly more damaged in all diagonals and rows. The value of TFE in the inner ring (IR) is 1.8 to 4.6 times higher than in the outer ring (OR), and similarly, the value of MDFE of the IR is 1.2 to 3.6 times higher than that of the OR. As described in the previous section, the values of the contact reaction and the contact angle are equal in magnitude for the inner and the outer ring contacts; therefore, no differences can be made in terms of load or motion, and consequently, the different damage can only be attributed to the different contact curvatures of the bearing raceway where the inner contact is convex with respect to the longitudinal curvature (r_b vs. r_pw) and concave with respect to the transversal direction (r_b vs. r_p), while the outer ring contact is concave with respect to both transversal and longitudinal curvatures. Therefore, the increasing main diameter of the bearings would aid in balancing the damage at the inner ring.

Figure 7 shows the density of the friction energy, ρ_F, along the different bearing raceways, i.e. the x axis, the longitudinal direction of the raceway, and the y axis, the transversal direction, according to the local coordinate system of Fig. 5. The plotted sectors correspond to the location of the maximum values at D2R1 from 225 to 320° denoted by the regions a and b for DLC 1.2 and for DLC 6.4 (see Fig. 6). In order to facilitate visualization and comparison between adjacent ball damage, the distance between marks was reduced. The results show notable differences in the damage shapes for every contact, and some look similar to previous experimental results (Cubillas et al., 2021b). All marks can be bounded in an elliptical area, which is notably more rounded at the OR and more flattened at the IR; this, as previously mentioned, is directly influenced by the contact curvatures of the rings. The forms of damage are significantly different for DCL 1.2 and DCL 6.4. In the case of DCL 1.2, there are two intensity zones within each footprint: a central, elliptical zone and a horseshoe-shaped zone at the rear of the tread. In contrast, the DCL 6.4 damage shows only one damage zone in the rear of the elliptically shaped tread. This fact may be related to the non-linear stiffness and how it affects the load ratio and contact angle. That is, for low loads, the contact is less stiff in relative terms as an increase in the load results in an increase in the contact angle. As the load increases, the stiffness of the contact is greater and the change in contact angle is smaller, and therefore the rolling is less. Under these conditions, the contribution from radial fretting damage, generated by the pulsating load, and from rotational fretting damage changes.

Figure 7 Density of dissipated friction energy, ρ_F, at regions a and b, according to Fig. 5.

3.3 On the effect of the wind speed

Figure 8 shows the results of the contribution of different wind speeds – 3, 5, and 7 m s⁻¹ for DLC 1.2 and 3, 11, and 25 m s⁻¹ for DLC 6.4 – to the total accumulated friction energy, TFE. This is unlike the previous analysis, where the sum of the contribution of each wind speed was presented. Due to the similarities of the damage distribution for both indicators, TFE and MDFE, only the results of TFE are presented. According to previous sections, the results are also classified by the load diagonal, row, and ring.

Figure 8 Results of the total dissipated friction energy for different wind speeds of the design load cases DLC 1.2 and DLC 6.4 for 10 min duration, where results are organized by load diagonal and row.

As mentioned, the results show the contribution of different wind speeds: 3, 5, and 7 m s⁻¹ for DLC 1.2 and 3, 11, and 25 m s⁻¹ for DLC 6.4. As can be observed, a similar distribution of the dissipated energy where the location of the damage seems not to be affected by the wind speed and the maximum values can be found at 60, 135, and 190° at D1 and at 30 and 270° at D2. Therefore, it can be concluded that the position of the critical zones is not dependent on wind speed. If the wind speed does not alter the location of the damage, it seems that the contribution of the fluctuating wind load is low. Therefore, it can be inferred that the main source of damage is caused by the swinging of the blade weight when the rotor is in motion.

While the location of the damage is invariant, a notable and expected effect of the wind speed is observed on the value of TFE, with the results showing a growing tendency of the dissipated energy with increasing wind speed. Figure 9 shows the normalized value of the maximum TFE and MDFE for the different wind speeds. Despite the analogous growing tendency for both indicators with the growing wind speed, the tendency is significantly different: while TFE has exponential growth, MDFE shows a logarithmic tendency. The TFE indicator refers to a global term after integrating local damage; i.e. it is a parameter that takes into account the total damaged area. In contrast, MDFE is a local parameter, and it takes into account the intensity of damage at a given point. Therefore, the different trends of the two damage parameters indicate different effects. MDFE indicates that the damage intensity grows logarithmically with wind speed, while the TFE parameter shows that the damaged area grows exponentially.

Figure 9 Evolution of the normalized values of the maximum density of friction energy (MDFE) and the total dissipated friction energy (TFE) as a function of the wind speed and the logarithmic and exponential regressions, respectively. (a) DLC 1.2; (b) DLC 6.4.

4 Conclusions

This work started with the aims of analysing the performance of fretting damage in productive and non-productive periods of pitch bearings under realistic time series and determining the critical locations and conditions that favoured the development of the damage. As a result, the following conclusions have been made:

• The energy-based wear model shows a critical region between 270 and 315° for the inner ring and at the first row, where the total accumulated friction energy (TFE) and the maximum density of friction energy (MDFE) are located in both operational and non-operational conditions.

• An analysis of the loads and contact angle distribution shows a noticeably better correlation between the friction energy and the variable load than the variable angle, and consequently, this implies that radial fretting is the main source of damage.

• The evaluation of the effects of the wind speed on the damage shows a similar distribution of damage where maximum values remain at the same level. Therefore, the critical locations do not vary with the intensity of the wind. However, important differences are found in the value of the accumulated energy indicators, showing exponential behaviour of TFE, and a logarithm trend of MDFE as a function of the wind speed. Therefore, the intensity of the damage is expected to be exponential with the wind speed.

• For time series of the same duration, the damage intensity obtained in productive times associated with the DLC 1.2 series is much higher than in the DLC 6.4 series. Consequently, the productive times are much more critically affected by damage development.

• The results obtained in this work shed light on the fretting damage mechanism in pitch bearings under different conditions of wind and operation. However, there is still a long way to go. Among these tasks, the most obvious would be to validate the proposed methodology on a real scale, with and without lubrication, to see the validation framework, as well as needing to implement additional formulations for modelling the effect of grease. Additionally, the formulation should be extended to the prediction of wear over time, making the methodology usable as a prediction tool for design.

Appendix A: Contact problem extended

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }{a}^{*}={\left({\displaystyle \frac{\mathrm{2}\cdot {\mathcal{H}}^{\mathrm{2}}\cdot \mathit{\xi }}{\mathit{\pi }}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(A2)}& {\displaystyle }& {\displaystyle }{b}^{*}={\left({\displaystyle \frac{\mathrm{2}\cdot {\mathit{\xi }}_i}{\mathit{\pi }\cdot {\mathcal{H}}_i}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(A3)}& {\displaystyle }& {\displaystyle }K={\displaystyle \frac{{\mathrm{2}}^{\mathrm{2}/\mathrm{3}}}{\mathrm{3}\cdot {\mathit{\delta }}^{*}\cdot \left(\mathrm{1}-v^{\mathrm{2}}\right)\cdot {\left(S_{\mathit{\rho }}\right)}^{\mathrm{2}/\mathrm{3}}}},\end{array}$$

where v is Poisson's ratio, S_ρ is the sum of the contact curvatures, and δ^* is the dimensionless contact deformation.

The sum of the contact curvatures is calculated as

$$\begin{array}{}\text{(A4)}& S_{\mathit{\rho }}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{w}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{w}}}}\pm {\displaystyle \frac{\cos \left(\mathit{\alpha }\right)}{r_{\mathrm{pw}}\pm r_{\mathrm{w}}\cdot \cos \left(\mathit{\alpha }\right)}}\pm {\displaystyle \frac{\mathrm{1}}{r_{\mathrm{p}}}},\end{array}$$

where r_pw is the main bearing radius, r_p is the raceway radius, α is the contact angle, and the sign of both terms is taken as positive for the outer contact and negative for the inner contact.

The dimensionless contact deformation is calculated as

$$\begin{array}{}\text{(A5)}& {\mathit{\delta }}^{*}={\left({\displaystyle \frac{\mathrm{4}\cdot F^{\mathrm{3}}}{{\mathit{\pi }}^{\mathrm{2}}\cdot k^{\mathrm{2}}\cdot \mathit{\xi }}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where ξ is the elliptical parameter of the first kind and F is the elliptical parameter of the second kind. These are calculated as

$$\begin{array}{}\text{(A6)}& {\displaystyle }& {\displaystyle }{\mathcal{H}}_i={\mathit{\rho }}^{\frac{\mathrm{2}}{\mathit{\pi }}},\text{(A7)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_i=\mathrm{1}+{\displaystyle \frac{\mathit{\pi }-\mathrm{1}}{\mathrm{2}\cdot \mathit{\rho }}},\text{and}\text{(A8)}& {\displaystyle }& {\displaystyle }F={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}+\left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}-\mathrm{1}\right)\cdot \ln \left(\mathit{\rho }\right),\end{array}$$

respectively, where ρ is the radius ratio, calculated as

$$\begin{array}{}\text{(A9)}& \mathit{\rho }={\displaystyle \frac{\mathrm{4}\cdot r_{\mathrm{pw}}\cdot c\cdot r_{\mathrm{w}}}{r_{\mathrm{w}}\cdot (\mathrm{2}\cdot c-\mathrm{1})\cdot \left(r_{\mathrm{pw}}\pm r_{\mathrm{w}}\cdot \cos \left(\mathit{\alpha }\right)\right)}},\end{array}$$

and c is the conformity aspect ratio, $c=\mathrm{2}{r}_{\mathrm{p}}/r_{\mathrm{b}}$; the sign is considered positive for the outer contact and negative for the inner contact.

Appendix B: Nomenclature

Abbreviations
IR	Inner ring
OR	Outer ring
FEM	Finite-element model
CoF	Coefficient of friction
TFE	Total accumulated friction energy
MDFE	Maximum density of friction energy
DLC	Design load case
DFE	Density of the friction energy

Symbols
a	Contact ellipse major semiaxis
a0	Distance from a point to the ellipse edge in
	the x direction
a*	Dimensionless major semiaxis
B	Bearing ring width
b	Contact ellipse minor semiaxis
b0	Distance from a point to the ellipse edge in
	the y direction
b*	Dimensionless minor semiaxis
C	Creepage coefficient
c	Conformity aspect ratio
Ds	Bolt diameter
E	Elastic modulus of the material
EF	Dissipated friction energy
ER	Rolling energy
F	Elliptical parameter of the second kind
G	Shear modulus
H	Bearing ring height
kw	Coefficient of wear
L	Flexibility parameter
Mx	Moment component x at bearing
My	Moment component y at bearing
P	Contact pressure
Q	Reaction force at the ball contact
rw	Radius of the ball
rpw	Bearing pitch radius
re,i	Cross-sectional groove radius, e for the outer
	and i for the inner ring
$r_{\mathrm{0}}^{\prime }$	Initial effective radius
r′	Effective radius
S	Local sliding
Sρ	Sum of the contact curvatures
T	Tangential stresses
t	Time
v	Poisson's ratio
Vw	Wear volume
W	Local rigid slip
WNL	Local rigid slip of the normal load
WRo	Local rigid slip of the rolling effects
x′, y′	Local coordinate system directions
z	Number of balls
ε	Rolling creepage
Tl	Traction limits of the contact
U	Local elastic slip
α	Contact angle
ρF	Density of friction energy
δ*	Dimensionless contact deformation
ℋ	Elliptical parameter of the second kind
ξ	Elliptical parameter of the first kind
μ	Value of the CoF
ΦX	Probability function of the normal distribution
X	Density of the normal probability distribution