Virtual sensing of wind turbine hub loads and drivetrain fatigue damage

This paper presents a Digital Twin for virtual sensing of wind turbine aerodynamic hub loads, as well as monitoring the accumulated fatigue damage and remaining useful life in drivetrain bearings based on measurements of the Supervisory Control and Data Acquisition (SCADA) and the drivetrain condition monitoring system (CMS). The aerodynamic load estimation is realized with data-driven regression models, while the estimation of local bearing loads and damage is conducted with physics-based, analytical models. Field measurements of the DOE 1.5 research turbine are used for model training and validation. The results show low errors of 6.4% and 1.1% in the predicted damage at the main and the generator side high-speed bearing respectively.


Introduction
Offshore wind turbine installations are projected to accelerate rapidly in the near future driven by better wind resources and higher social acceptance rates compared to onshore sites [24]. However, a major economic limitation of off- shore wind turbines are high operational and maintenance expenditures (OPEX), which amount to about 34 % of the levelized cost of energy (LCOE) [19]. These are caused by lower reliability due to harsher environmental conditions and time-consuming replacement or repair due to difficulties accessing the site and dependency on good weather conditions. A major contributor to OPEX is the geared drivetrain with frequent failures and long downtimes and is thus the subject of current research [23].
Digital Twin (DT) is an emerging technology with prospects of decreasing operational and maintenance expenditures and improving the market competitiveness of offshore wind farms. The wind turbine drivetrain DT proposed by the authors in [readacted] would enable monitoring of fatigue loads at otherwise inaccessible locations such as bearings and gear contacts using 'virtual sensors' and thus support Remaining Useful Life (RUL) assessment based on the true load history. A framework with the three components Data, Virtual model and Decision support is envisioned for this objective (Fig. 1). Data is collected continuously by sensors of the Supervisory Control and Data Acquisition system (SCADA), the condition monitoring system (CMS), and from other sources such as metocean forecasts and maintenance logs. The Virtual model comprises of decoupled simulation models to represent the physics at different scales. Aerodynamics and structural dynamics of tower and blades are captured with the aeroelastic model. The resulting hub loads and nacelle motions are imposed as boundary conditions on the drivetrain model, which produces local bearing and gear forces. These are used as input for the fatigue damage model consisting of stress cycle counting algorithms and S-N curves. The virtual model and its physical counterpart are synchronized with real-time field measurements using state estimators such as the Kalman Filter. The synchronization, also referred to as data fusion or Digital Twinning, is essential as it facilitates measurements of virtual sensors in the synchronized model. The virtual sensor measurements are converted to value adding information for the turbine operator in the component called Decision support. In this study the focus is on RUL assessment of drivetrain components, which is necessary to move from corrective to predictive maintenance strategies.
Preliminary investigations on the feasibility of the proposed DT have been conducted in a numerical case study in [redacted]. In this study the proposed DT is further validated with with field measurements of the DOE 1.5 MW turbine instrumented by NREL [18]. Main bearing loads estimated with the proposed virtual sensing method are compared to field measurements obtained from shaft strain gauges under different operational conditions.
Other studies are found in literature, which are concerned with estimating main shaft loads with a virtual sensing approach. Several works pursue an inverse approach for real-time estimation of the rotor torque based on SCADA measurements. By simplifying the drivetrain dynamics to a two degrees of freedom (DOF) torsional system, the SCADA signals generator torque along with the LSS and HSS speeds contain enough information to predict the unknown rotor torque [3,12,15,17]. The drivetrain model may be constructed without knowledge of manufacturer's specifications using data-driven system identification techniques such as least-squares estimators [15]. Data fusion is realized with state estimators, for example Kalman Filters (KF) [15], or regularization methods [17].
Several works are concerned with virtual sensing of rotor thrust, but with the objective of structural health monitoring of the tower rather than drivetrain components. A common approach involves constructing a dynamic, flexible tower model and use state estimation methods based on tower top acceleration and/or strain gauge measurements [20].
Notable publications that fit the proposed Digital Twin framework ( Fig. 1) are presented by Branlard et al. and Azzam et al. [2,4]. Branlard et al. present a holistic wind turbine DT capable of estimating both thrust and torque based on SCADA measurements and a linearized aeroelastic model [4]. Validation of the DT is conducted with both simulation and field measurements, however the scope is limited to structural dynamics of the tower and blades. Azzam et al. present a DT that also considers drivetrain dynamics [2]. The DT is constructed by regression on aeroelastic and drivetrain multi-body simulations and serves the purpose of estimating all six main shaft load components based on SCADA measurements. Unfortunately, their work is purely numerical and not supported by experimental or field measurements.
The novel contributions of this paper in comparison to existing research are summarized as follows.
Validating the concept of virtual sensing of drivetrain loads with field measurements rather than numerical simulations Monitoring fatigue damage in drivetrain components rather than the tower or blades Leveraging high-frequency CMS vibration measurements as opposed to using only SCADA measurements The remainder of this paper is organized as follows. Sect. 2 presents the proposed methodology for virtual sensing of drivetrain loads and remaining useful life estimation. It follows a discussion of the errors between field measurements and predictions using the virtual sensors in Sect. 3. Concluding remarks are given in Sect. 4.

Methodology
The overall methodology is illustrated in Fig. 2

Description of dataset
The dataset was acquired as part of a field measurement campaign with the U.S. Department of Energy 1.5 MW (DOE 1.5) turbine at the National Renewable Energy Laboratory (NREL) [18]. The DOE 1.5 turbine is based on a commercial GE 1.5 SLE turbine with a custom configuration. In addition to a standard SCADA system and drivetrain CMS the turbine is equipped with strain gauges at the tower base, tower top, blade roots and the main shaft to fully monitor multiaxial aerodynamic loads. The sensor signals used in this study are listed in Table 1. The following SCADA signals are considered in this study, which are reportedly sensitive to the main shaft loading: Active power, HSS and LSS speed, Nacelle wind speed, as well as tower top acceleration. The CMS sensors are installed in a typical configuration and positioned on the housing of the generator (Gen), the high-speed gear stage (HSS), planetary gear stage (PL) and each of the torque arms (TA). The aerodynamic loads at the rotor hub including the torque, pitch moment, yaw moment and thrust are  measured with strain gauges at the main shaft downwind of the main bearing and at the tower base. The calibration of the strain gauges is described in [18]. The dataset used in this study comprises of a total of 830 measurements of 10 min length recorded from 31. Oct 2018 to 05. Dec 2018. The sampling frequency is 50 Hz for all signals.

Data postprocessing
The dataset is filtered for normal power production, which is identified by three criteria Main shaft speed > 10.5 rpm Blade 1 pitch angle < 50 ı Active power > 0 kW In addition, a moving average filter with window size of 1 s is applied on the recorded strain gauge signals.
Best practice in drivetrain condition monitoring is the extraction of statistical features, which are indicative of faults or damage [16]. The recorded SCADA and CMS measurements are partitioned into 10 min segments and the features listed in Table 2 are calculated for each segment. These include a wide range of the most commonly used features in the time domain (x) and frequency domain (S xx ). The features that are eventually utilized as input for the regression models are determined by a sensitivity analysis in Sect. 3.2.

Data-driven regression models
Regression models are used in this study to map the predictors, the SCADA and CMS statistical features, onto the targets, the aerodynamic hub loads. Several linear and nonlinear regression models are investigated for this purpose including Linear Regression (LR), Support Vector Machine (SVM) and Tree ensembles, as described in Table 3. For a detailed description of each regression model type it is referred to [6]. Implementation and training is realized with MATLAB's Statistics and Machine Learning Toolbox [9]. The dataset is partitioned 80/20 into training data and test data, and the models are regressed onto the training data

Physics-based drivetrain model
The DOE 1.5 MW turbine is instrumented with strain gauges at the blade roots, the main shaft and the tower top and base to monitor the multiaxial loading of the turbine. The aerodynamic loads at the rotor hub, as well as loads at the main bearing and HSS-GS bearing are calculated from strain gauge measurements using an analytical model, presented in [1,5]. The analytical model assumes steady state operation and neglects any torsional or bending dynamics of the main shaft and the tower. With this assumption it is possible to determine bearing loads by moment balances.
First, the aerodynamic moments including torque M a,x , pitch moment M a,y and yaw moment M a,z are determined from the measured main shaft moments M ms by moment balance around the main bearing (Fig. 3) and expressed in the fixed coordinate frame at the hub. The thrust F a,x is calculated from the tower base moment M tb,y by moment balance around the tower base (Fig. 3) where M ms;y0 and M tb;y0 are gravitational moments due to the rotor overhang expressed at the main bearing and tower base respectively, h denotes the tower height, d 1 the distance from the hub to the tower top and˛the main shaft tilt angle (Table 4) The installed main bearing is a SKF 240/600 BC spherical roller bearing in a 3-point configuration and thus considered to only support radial and axial forces. The torque arms are also considered to only experience radial and axial forces and the stiffness of the generator coupling is neglected. In steady state operation the main bearing forces F mb are then calculated as follows (Eqs. 5-7). The axial force is governed by thrust, while the radial force is governed by the yaw and pitch moments.
where F mb;x0 and F mb;z0 is the rotor, shaft and gearbox weight projected onto the x an z axis respectively and d GB is the distance from the main bearing to the torque arms ( Table 4). The HSS-GS bearing is a SKF NU232 cylindrical roller bearing and thus only supports radial forces. The radial force is governed by the transmitted gear force at the highspeed gear stage, which is calculated from the rotor torque by neglecting all torsional dynamics where i GB denotes the gearbox ratio, r b the base radius of the pinion and d RS , d GS the distance from the generatorand rotor side HSS bearings to the pinion center (Table 4).

Fatigue damage and remaining useful life
The bearing fatigue damage and remaining useful life is based on ISO 281 [7], which defines the equivalent dynamic load P for cylindrical roller bearings (CRB) and tapered roller bearings (TRB) as for CRB: P = F rad (10) for TRB: where Y 1 , Y 2 , e are bearing specific parameters ( Table 4). The equivalent dynamic load is calculated with 10 min average load estimates denoted as P i . For each 10 min section indexed by i the permissible stress cycles N i are then calculated with the bearing lifetime equation where C is the basic dynamic load rating and m equals 10/3 for roller bearings. The experienced stress cycles n i are determined using the load duration distribution (LDD) method, which counts one stress cycle per shaft revolution due to entering and exiting the bearing load zone [13]. Using 10 min average shaft speeds ! i the LDD method simplifies to where t equals 10 min. It follows the dimensionless shortterm fatigue damage D ST i , which is defined as the ratio of experienced to permissible stress cycles The long-term damage D LT .t/ is obtained with the Palmgren-Miner linear damage hypothesis by summation of all previous short-term damage and is updated in 10 min intervals for real-time damage monitoring By definition, the bearing has consumed its damage reserves and reached its end of life at D LT = 1. With a nominal life t nom of 20 years the remaining useful life RUL is then calculated as follows

Damage uncertainty model
Using 10 min average load estimates for the damage calculation reduces computational costs and enables real-time monitoring, however it introduces uncertainties by neglecting high-frequency load fluctuations, which may originate in the aerodynamics or internal drivetrain dynamics. The damage is generally underestimated with averaged loads, since load peaks are disproportionally more damaging than load minima due to the exponentiation with m (Eq. 12). The uncertainty avg is expressed as the ratio of the true shortterm damage D ST 50 Hz measured at 50 Hz and the short term damage based on 10 min average load estimates D ST 10 min-avg . The fluctuations of the equivalent dynamic load within a 10 min period are modelled with a statistical variable X N. ; /, which is normally distributed with mean value and standard deviation . It is further assumed that variations of the shaft speed are negligible, such that Eq. 13 remains valid. It follows for the uncertainty avg where the expected values E are given by the law of the unconscious statistician (LOTUS) [8] using the standard It is evident that the uncertainty avg is only a function of the 10 min mean and standard deviation, which are both estimated with data-driven regression models (Sect. 2.3).

Measured hub loads and fatigue damage
Shown in Fig. 4 are the distributions of the measured hub loads for verification of the results. The mean torque follows the analytical thrust curve and levels out at rated torque, while the highest torque variance is observed at about 10 m/s slightly below rated wind speed of 14 m/s. This behaviour is similarly reported in other works [13] and is likely due to effects of the pitch controller, which frequently activates and deactivates in this region causing high torque amplitudes. The aerodynamic pitch moment is predominantly a result of thrust differences between the upper and lower rotor disk due to the vertical wind profile. Positive trends of the mean and variance with reference to wind speed is observed. The yaw moment is similarly a result of aerodynamic imbalance, predominantly yaw misalignment. Contrary to the pitch moment, the yaw moment is centered around zero mean and is independent of wind speed. The measured thrust agrees well with simulated thrust curves, as demonstrated in [5]. The highest variance in thrust is measured at around 8 m/s, which is slightly lower than the peak of torque variance. The calculated bearing damage based on the measured hub loads are presented in Fig. 5. It is emphasized here that rotating machine elements such as bearings and gears experience stress cycles even at stationary environmental loads due to the shaft rotation. For this reason the LDD method [13] is used in this study for stress cycle counting as opposed to the rainflow counting method commonly used for (non-rotating) structural elements. The fatigue damage at the main bearing comprises of two components, an axial component D ax = XF ax =P D due to thrust and a radial component D rad = Y F rad =P D due to gravity and pitch moments. The maximum in fatigue damage is observed at 11 m/s and coincides with the thrust peak. In this operational region the axial forces dominate and amount to about 66% of the equivalent dynamic load P . At wind speeds above 16 m/s the contribution of radial forces due to pitch moments becomes dominant and below 8 m/s with relatively low aerodynamic loads the contribution of gravity forces becomes dominant.
The HSS-GS bearing experiences only radial forces, which are considered to be proportional to the rotor torque (Eq. 8). Thus, the fatigue damage is governed by the mean rotor torque and reaches its maximum above rated wind speeds.

Sensitivity analysis
A sensitivity analysis is conducted with the objective of dimensionality reduction of the predictor variables. The sensor signals (Table 1) and statistical features (Table 2) are selected, which are the best predictors of hub loads according to the metric of the Neyman-Pearson correlation coefficient. Presented in Fig. 6 are the ten best performing signals SCADA signals (in blue) and CMS signals (in red) for each hub load component.
The generator torque is as expected an excellent predictor of both the mean and the standard deviation of rotor torque with R > 0.99. Prediction of the absolute values of the bending moments on the other hand is challenging, as neither of the SCADA or CMS signals show statistically significant correlation (R < 0.5). However, the moment standard deviations show high correlation (R = 0.88) with respect to the wind speed, as well as other SCADA signals. The torque correlates well with all SCADA signals, as well K Fig. 7 Normalized RMSE between measured and predicted hub loads with different regression models and different SCADA/CMS input as the CMS vibrations at the HSS, the generator and the nacelle frame.

Hub load estimation
Several regression models, as described in Table 4, are used to map sensor measurements onto the aerodynamic hub loads. Two different scenarios of sensor input are considered: (1) only SCADA signals, (2) combined SCADA and CMS signals. This serves the purpose of assessing the added value of CMS vibration data and validating the novel approach of virtual sensing based on vibration measurements. The metric for model performance is the root mean square error (RMSE) between measured and predicted loads using 5-fold cross validation. Shown in Fig. 7 is the RMSE normalized to the maximum value of each hub load.
It is evident that the estimation of the mean and standard deviation of the rotor torque is accurate with minimum RMSE of 0.24% and 0.46% respectively. The best performance is observed is observed with a simple linear Concerning the bending moments, it appears that estimating the mean value is much more challenging than estimating the standard deviation. A possible reason is that the mean bending moments unlike torque and thrust do not show a clear trend with respect to wind speed (Fig. 4). The inclusion of CMS vibration data slightly improves the prediction accuracy of bending moments in most cases. Nonlinear regression models are preferable, since the relationship of bending moments and dynamic drivetrain responses appear to be non-linear.
The mean thrust as well as the standard deviation is estimated with relatively low RMSE of 8.0% and 6.0%. It is clear that non-linear regression models are necessary to capture the non-linear behaviour such as the thrust-wind speed curve (Fig. 4). In this case CMS vibration data do not appear to increase performance.

Fatigue damage and remaining useful life
The measured and estimated hub loads discussed in the previous section are converted into short-term (10 min) fatigue damage in the main bearing and the HSS-GS bearing using Eqs. 5-16. The RMSE of the fatigue damage normalized to its maximum value is presented in Fig. 8.
The damage in the main bearing is estimated with high accuracy (RMSE = 6.4%) despite high uncertainty in estimating the bending moments. These results suggest that the damage in the main bearing is governed by thrust, which can be estimated more accurately. The best performance is achieved by the quadratic SVM, which is able to capture the non-linear behaviour best. For monitoring the damage in the HSS-GS bearing a linear regression model suffices, which results in RMSE of 1.1%.
It appears that the inclusion of high-frequency CMS vibration measurements does not provide much value for monitoring bearing fatigue damage and that the considered 10 min average SCADA measurements are sufficient to estimate the damage within a 6.4% error margin. Fig. 9 presents the measured and predicted RUL with the best performing regression model. During the recorded time frame of 138.3 h the measured RUL of the main and HSS-GS bearing is reduced only by 20.8 h and 17.5 h respectively. The discrepancy can be attributed to conservative design, for example in the selection of design load cases (DLC), which are more severe than the actual experienced environmental conditions. Furthermore, the sample size is relatively small and the time frame of the recordings of 31. Oct to 05. Dec is not representative for seasonality effects.
The RUL is overestimated significantly despite high accuracy in the predicted 10 min average loads. This is caused by high-frequency load dynamics for example from turbulence or internal drivetrain dynamics, which are not accounted for with 10 min average load estimates. The discrepancy is partially compensated with the damage uncer-tainty model (Sect. 2.6). A good agreement with the measured RUL is observed at the HSS-GS bearing, while at the main bearing there remains a larger error possibly due to higher uncertainties in predicting bending moments and thrust.

Method limitations and further work
While the presented Digital Twin exhibits high accuracy in predicting aerodynamic loads and bearing damage, it is crucial to discuss the method assumptions and associated uncertainties, which limit the applicability of this method.
Field measurements: The data used in this study (Sect. 2.1) are representative for commercial SCADA and CMS measurements with the exception of the wind speed. The wind speed data were acquired with a MET mast about 150 m downwind of the turbine. Commercial wind turbine SCADA systems, however, mostly rely on nacelle mounted anemometers, which suffer from greater inaccuracies due to wake effects. The additional measurement uncertainty can be estimated with a coefficient of variation (COV) of 1 − 3% according to Toft et al. [21].
Aeroelastic model: The presented regression model (Sect. 2.3) relies on a training data set of aerodynamic loads, in this case obtained by strain gauge measurements, which are not available in commercial wind turbines. Alternatively, it is possible to emulate field measurements with measurements from high-fidelity simulation models, similar to the approach of Azzam et al. [2]. Naturally, this shifts the challenge to the model construction and validation and introduces additional uncertainties due to modelling errors. Such uncertainties can be approximated with a COV of 5% according to Nejad et al. [13], however it is difficult to make generalized statements. In further studies it is planned to compare the data-driven regression models with an aeroelastic model of the DOE 1.5 turbine, which has been developed and validated by other authors [5], in order to quantify modelling uncertainties.
Drivetrain model: State-of-the-art drivetrain models are highly complex multibody simulation (MBS) models [14,22], and are not suitable as Digital Twin models, as discussed in [11]. First, the high number of degrees of freedom (DOF) make them numerically expensive and not capable of real-time simulation, which is necessary for online monitoring purposes. Secondly, wind turbine operators do not have the means of developing and validating complex models, since the drivetrain specifications are largely confidential to the OEMs. For this reason, a relatively simple drivetrain model is used in this article, which assumes a quasi-static transmission of torque and neglects all internal dynamics including effects of component flexibility, multi-body interaction and excitations from gear meshing or roller bearings (Sect. 2.4). The effects of internal dynamics on bearing fatigue damage are expected to be relatively small, as suggested by the results of a previous numerical case study [10], where RMSE of 5-15% in the bearing fatigue damage were observed. However, the scope of the numerical case study was limited to the high-speed gear stage bearings, normal power production at rated wind speed and one drivetrain configuration. Further numerical investigations are scheduled better quantify the modelling errors of such reduced order drivetain models.

Conclusion
This paper presents a Digital Twin for virtual sensing of wind turbine hub loads based on SCADA and CMS measurements, as well as monitoring the accumulated fatigue damage and remaining useful life in the main and HSS-GS bearing. The Digital Twin is constructed for the DOE 1.5 research turbine [18] and evaluated with field measurements. Several data-driven regression models including linear regression models, support vector machines and tree ensembles are trained on field measurements for the aerodynamic hub load estimation. For calculation of local bearing loads a low-fidelity physics-based model is constructed with the assumption of steady-state operation. The remaining useful life is calculated based on the consumed fatigue damage reserves according to ISO 281 [7].
While the estimation of rotor torque and thrust is accurate with RMSE of 0.24% and 6.0%, it proves to be much more challenging to estimate the yaw and pitch bending moments. The measured bending moments appear to be highly stochastic and do not show statistically significant correlation (R < 0.5) with any of the considered SCADA and CMS measurements.
Nonetheless, relatively low RMSE of 6.4% in the 10 min fatigue damage are observed at the main bearing despite the high uncertainty in the bending moment estimates. It appears that the damage in the main bearing is governed by thrust, which is estimated much more accurately than the bending moments. The damage at the HSS-GS bearings is assumed to only depend on the drivetrain torque and can thus be estimated with high accuracy (RMSE = 1.1%).
The main contribution of this article is the knowledge transfer of the virtual sensing concept from wind turbine structural elements to drivetrain components, and validation of the concept with field measurements. With regards to the quality and availability of physical sensor measurements the proposed virtual sensors are feasible. SCADA and CMS data contain sufficient information for accurate monitoring of bearing fatigue damage. Challenges are identified in the multi-body drivetrain dynamics, which are much more complex than the dynamics of the tower and blades. However, developing and validating models to capture complex drivetrain dynamics is difficult based on the information that is available to wind turbine operators. Low fidelity, quasistatic models, which largely neglect internal drivetrain dynamics, are shown to produce low errors of bearing fatigue damage, and are thus proposed for virtual sensing purposes. Further investigations are planned to quantify the uncertainties introduced by quasi-static drivetrain models.
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital).
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/.