## Abstract

In this article a virtual sensor for online load monitoring and subsequent remaining useful life (RUL) assessment of wind turbine gearbox bearings is presented. Utilizing a Digital Twin framework the virtual sensor combines data from readily available sensors of the condition monitoring (CMS) and supervisory control and data acquisition (SCADA) system with a physics-based gearbox model. Different state estimation methods including Kalman filter, Least-square estimator, and a quasi-static approach are employed for load estimation. For RUL assessment the accumulated fatigue damage is calculated with the Palmgren–Miner model. A case study using simulation measurements from a high-fidelity gearbox model is conducted to evaluate the proposed method. Estimated loads at the considered intermediate and high-speed shaft bearings show moderate to high correlation (R = 0.50 − 0.96) to measurements, as lower frequency internal dynamics are not fully captured. The estimated fatigue damage differs by 5–15% from measurements.

## 1 Introduction

Recent market trends show an increased shift toward offshore wind turbine installations due to the higher energy yield and fewer issues with land displacement and noise [1]. However, offshore sites face additional reliability challenges. Replacement or repair of components is expensive and time consuming due to difficulties accessing the site and dependency on good weather conditions. Thus, unscheduled down times as a result of component failure can lead to high operational and maintenance expenditures (O&M). For offshore wind turbines the O&M expenditures can reach 34% of the levelized cost of energy (LCOE), twice as much as for land-based turbines [2]. A major contributor to the O&M expenditures is the gearbox with a failure rate of 0.1–0.15/year and average downtimes of 6 days per failure [3,4].

*Predictive maintenance* strategies are proposed in the offshore industry to increase reliability and availability, and decrease O&M expenditures. As a subcategory of condition-based maintenance (CBM) predictive maintenance depends on continuous monitoring of the systems’ operational condition for the assessment of the remaining useful life (RUL). Alerts are triggered in the case of severe deviation of RUL to nominal life and the operator may schedule immediate maintenance tasks in addition to regular, time-based maintenance routines. Currently, the predictive capabilities of condition monitoring systems (CMS) are limited. In practice, trends of vibration-based statistical features are extrapolated to a predefined failure threshold to predict the RUL, which suffers from a high level of uncertainty [5].

*Digital twin* (DT) is identified as an emerging technology that could facilitate predictive maintenance strategies [6]. DT can be described as a virtual representation of a physical asset enabled through real-time measurements and simulators for the purpose of improved decision making [7]. The authors previously proposed a Digital Twin framework, shown in Fig. 1, with the three components *Virtual model, Data*, and *Decision support* to move toward predictive maintenance [8].

*Virtual models* of wind turbines have matured in the past two decades to a high level of fidelity. Generally a decoupled approach is employed with aeroelastic models for global dynamics, multi-body simulation (MBS) models for drivetrain dynamics, and finite element (FE) models for component dynamics [9]. While many authors associate DT with high-fidelity models, recent publications move toward reduced-order models (ROM) to meet requirements of computational speed for real-time monitoring tasks [8]. The DT model is updated with data such that it virtually experiences the same environment as its physical counterpart.

*Data* that can be leveraged in wind turbines are sensor measurements of the drivetrain CMS or the supervisory control and data acquisition (SCADA) system. Typical signals include vibration on the gearbox housing, electrical signatures of the generator and shaft speeds.

*Decision support* is a collective term for services that the DT provides to assist the operators’ maintenance or control decisions. Focus of this research are methods for online monitoring of loads in drivetrain components (gears, bearings) and subsequent RUL estimation. Direct measurements of component level loads are difficult and require custom solutions, such as bearings with integrated strain gauges, which are generally not available for commercial wind turbines. Hence, indirect (or inverse) methods of load estimation that combine more accessible sensor measurements and a DT model would provide a cost-effective alternative. This procedure is often referred to as *Virtual sensing*, as it can be interpreted as taking measurements in a fully synchronized virtual model. Synchronization is achieved by continuously estimating the dynamic states of the system. Different state estimation methods are employed for this purpose, most prominently the Kalman filter and its variations, as well as least-square estimators.

The DT and virtual sensing approach are often pursued for estimating damage equivalent loads or stresses for structural health monitoring (SHM), for example in wind turbine towers. Virtual sensing of tower loads based on a limited number accelerometers and strain gauges has been demonstrated both in numerical and experimental studies [10,11]. However, limited research has been conducted for drivetrains with the exception of a study by Bosmans et al., who present a virtual sensor for wind turbine planetary gear loads based on strain gauge measurements and a FE model [12]. The use case of drivetrains comes with unique challenges: The internal dynamics of drivetrains are much more complex due to multi-body interactions and there are limitations in the existing drivetrain sensors (SCADA, CMS) related to signal resolution, sensor locations, and noise that make it difficult to observe the current dynamic state.

The main contribution of this work is to apply a virtual sensing approach that has proven to be effective in other areas to the use case of wind turbine drivetrains and demonstrate the proof of concept in a numerical case study. In addition, this article discusses the challenges that are unique to drivetrains and provides some recommendations on suitable sensor signals and state estimation methods.

The remainder of this paper is organized as follows. In Sec. 2, the mathematical development of the virtual load sensor and its use case in a wind turbine high-speed gear stage are outlined. Section 3 discusses the performance of the virtual sensor in a numerical case study. Lastly, Sec. 4 provides some concluding remarks.

## 2 Methodology

The virtual load sensor is developed in a DT framework with the three components of Data, Virtual model, and Decision support (Fig. 1). Section 2.1 presents the high-fidelity model of a reference wind turbine gearbox, which is linearized for integration with the virtual sensor, as shown in Sec. 2.2. Synthetic CMS and SCADA data are generated by means of simulation with high-fidelity models (Sec. 2.3). Different state estimators including the Kalman filter, Least-squares estimator, and a quasi-static approach are used for virtual sensing of bearing loads (Sec. 2.4). Subsequently, the accumulation of fatigue damage is tracked with the standard Palmgren–Miner model and the bearing lifetime equation according to ISO 281 (Sec. 2.5).

### 2.1 High-Fidelity Models.

A reference gearbox based on the NREL offshore 5 MW baseline wind turbine and mounted on the floating OC3 Hywind spar structure is considered in this study [13,14]. The reference gearbox was developed by Nejad et al. with reference to minimal weight and following offshore wind turbine design codes [15]. The gearbox comprises of two planetary and one parallel gear stage totalling to a gear ratio of 1:96.354. The main shaft support is a 4-point design with two main bearings to minimize non-torque loads entering the gearbox. A decoupled approach is employed, which is best practice for drivetrain simulation [9]. The global response to a set of environmental conditions is determined with the global model, which is implemented in the aero-hydro-servo-elastic code SIMO-Riflex-AeroDyn. The internal dynamics are then simulated with a high-fidelity gearbox model implemented in the multi-body simulation environment SIMPACK. External loads (torque and non-torque) are applied on the main shaft, the nacelle movements are applied on the bed plate and the generator torque is applied on the high-speed shaft to control the generator speed.

### 2.2 Linearized Model.

**M**denotes the diagonal mass matrix comprised of inertia terms. The stiffness and damping matrices contain terms from the elastic couplings of bearings (

**C**) and gear meshing (

_{b}, K_{b}**C**). The detailed matrix composition is given in Ref. [16]. External forces and moments crossing the system boundary at the generator and rotor side shaft interfaces are represented by

_{m}, K_{m}**f**

_{ex}∈

*R*

^{12×1}. The equations of motion are first linearized and then transformed into a set of first-order differential equations, the so-called state-space representation (Eq. (2)). In this step the time-variant mesh stiffness is reduced to a constant value $Cm\xaf$, hence the linearized model is unable to reproduce periodic excitation at the gear meshing frequency.

**x**of the state-space model is a stack of body-fixed displacements and velocities (Eq. (3)), while the external forces

**f**

_{ex}are split into known input variables

**u**(Eq. (4)) and unknown disturbance forces regarded as process noise

**w**(Eq. (5)). Of the 12 external force terms only the generator torque is considered available from SCADA measurements and thus a known input variable, while the remaining non-torque loads are modeled as white gaussian noise with covariance

**Q**. The system matrix

**A**(Eq. (6)) describes the dynamic state-transition and is obtained by rearranging mass, stiffness, and damping matrices [17]. The control matrix

**B**(Eq. (7)) represents the influence of input variables on the dynamic states.

**y**(Eq. (9)) are measurements of virtual accelerometers placed at the shaft bearings in axial and radial direction (Fig. 2) that represent CMS vibration sensors. These are analogously related to the state and input variables through a linear model (Eq. (8)), where

**C**denotes the observation matrix and

**D**the feedthrough matrix. The exact matrix composition cannot be shown, as these are generated numerically by SIMPACK’s linearization solvers.

**v**(Eq. (10)), which is modeled as white gaussian noise with covariance

**R**. In this case study measurement noise is added to the (exact) simulation measurements in postprocessing

**f**the general state-space model is augmented with an algebraic equation, which relates bearing loads to system states with the matrix

**K**(Eq. (11)). Since the bearings are considered as spring–damper elements in the drivetrain model, this relationship is linear. The matrix

**K**contains terms of bearing stiffness

**K**and damping

_{b}**C**and is generated numerically by SIMPACK.

_{b}**C, D, K**of the discrete model remain unchanged, as they only appear in algebraic equations, whereas

**A**can be derived as follows [17]:

_{d}, B_{d}**A, B, C, D, K**are calculated with SIMPACK’s built-in linearization solvers and integrated in the virtual load sensor in matlab.

### 2.3 Sensor Data.

In this numerical study simulation measurements from high-fidelity models are used to evaluate the proposed load estimation method. A reference load case at rated wind speed of 12 m/s (load case EC4, spar in [9]) is selected, since conditions near rated wind speeds are shown to induce the most severe drivetrain loads and have the highest contribution to long-term fatigue damage [18]. Six simulations each with a duration of 3800 s are conducted to comply with IEC 61400 guidelines [19]. The first 200 s are disregarded to avoid simulation start-up effects and the simulation time-step is set to 1 ms to capture high-frequency gear meshing dynamics. From the simulation results the generator torque, shaft vibration, and bearing loads are of interest. The generator torque and the shaft vibration are used to generate synthetic SCADA and CMS data as input for the load estimation method. Vibration signals are measured by virtual acceleration sensors mounted on the intermediate (IMS) and high-speed shaft (HSS) with a sampling frequency of 1 kHz. To capture yaw and pitch movements each shaft is equipped with two virtual sensors measuring axial and radial acceleration. White gaussian measurement noise $v\u223cN(0,R)$ is added to all acceleration measurements in postprocessing, where the covariance **R** is chosen, so that the signal-to-noise-ratio (SNR) is equal to 10 for all measurement signals. Additionally, the radial and axial loads at the IMS and HSS-bearings are extracted from simulation measurements for comparison with the estimated loads (see Table 1).

### 2.4 Virtual Load Sensor.

*n*. Three different state estimators with different levels of fidelity and requirements to measurement inputs are studied, Kalman filter, Least Squares, and Quasi-static. In each case, the system states $x^n$ are estimated first. Subsequently, the bearing loads are determined, as these are linearly dependent on the system states

#### 2.4.1 Kalman Filter.

**u**. The disturbance forces and moments

_{n}**w**on the system are not included in the prediction step, as they are regarded as process noise. The a priori estimated covariance $P^n\u22121|n\u22121$ is also predicted based on previous knowledge and the known process noise covariance

**Q**

**y**resulting in the a posteriori state estimates

_{n}**x**

_{n|n}**M**, which relates the confidence in state predictions of the physical model to the confidence in the measurement. With a high confidence in the physical model ($P^n|n\u22121\u21920$) the Kalman gain approaches zero; hence, the measurements update is given a low weight. On the other hand, with a high confidence in the measurements (

_{n}**R → 0**) the Kalman gain approaches

**C**. In this case the measurements have a higher significance compared to state predictions.

^{−1}#### 2.4.2 Quasi-Static Approach.

#### 2.4.3 Least-Squares Approach.

**C**

^{+}### 2.5 Fatigue Damage Model.

*n*

_{i}denotes the experienced stress cycles,

*N*

_{i}the number of cycles until failure, and

*i*indicates the stress range

*N*

_{i}the nominal bearing life equation according to ISO 281 [28] with the basic dynamic load rating

*C*and the equivalent bearing load

*P*is used

*P*is a linear combination of the axial and radial load with the factors

*X*and

*Y*, which are bearing specific values taken from the manufacturer’s data

*n*

_{i}are counted with the load duration distribution (LDD) method according to IEC 61400-4 [19]. The LDD method is applicable for rotating machinery components under slowly varying loads that experience cyclic loading due to entering and exiting the load zone each rotation. One stress cycle is counted for each rotation with a stress range equal to the current load.

*D*(

*t*) over time and observing the damage reserves. By definition, the end of the component’s nominal design life

*t*

_{nom}is reached at

*D*= 1 [18]. Hence, the RUL is retrieved as follows:

## 3 Discussion

The virtual load sensor using the state estimators Kalman filter (KF), Least-Square (LS), and Quasi-Static (QS) presented in Sec. 2.4 are evaluated in a case study. Radial loads at the IMS and HSS bearings estimated with the virtual sensor are compared to simulation measurements obtained from the high-fidelity drivetrain model outlined in Sec. 2.3. First, the correlation of estimated and measured loads is analyzed in the time and frequency domain. Secondly, the error in calculated fatigue damage is discussed.

### 3.1 Estimated Loads.

For a qualitative assessment the time series of measured and estimated radial loads in the IMS-A and HSS-A bearings is shown in Figs. 3 and 4. The measured loads *f* can be characterized as highly dynamic with lower frequency dynamic components (<10 Hz) as a result of slowly changing environmental conditions, as well as higher frequency dynamics (>100 Hz) induced by internal gearbox excitations such as gear meshing. The load estimates of the QS method $f^QS$ are sufficient to capture the long-term trend of bearing loads, but unable to reproduce any high-frequency internal dynamics seen at a time scale of 1 s. However, a slight bias at the HSS-A is observed, which could potentially be due to non-torque loads at the high-speed gear stage, which the torque-proportional QS does not take into account. The LS method produces load estimates $f^LS$ with high-frequency oscillations of similar amplitudes to measured loads, however at the IMS-A there appear to be several outliers, which significantly overestimate measured loads. In the low-frequency range, the LS method is not able to fully capture the internal gearbox dynamics. This is especially noticeable at the IMS-A, where the measured loads have a high-energy frequency component of about 5 Hz. The KF load estimates $f^KF$ are smoother and do not suffer from extreme outliers. Similar to the LS method high-frequency oscillations are captured well, while some lower frequency components are not reflected.

For the analysis of the behavior in the frequency domain the power spectral densities (PSD) of measured and estimated bearing loads are calculated, as shown in Figs. 5 and 6. The measured load spectrum shows several lower-frequency peaks (<10 Hz) and higher-frequency peaks at 80 Hz for IMS-A and at 464.2 Hz for HSS-A. The higher frequency peaks coincide with the gear meshing frequencies of the parallel and second planetary gear stage, respectively. The lower frequency peaks are not fully identified as of now.

The QS method matches the measured load spectrum of HSS-A reasonably well with the exception of the high-frequency range with the gear meshing peak, which is underestimated significantly. In the low frequency range the peaks at 4.75 Hz, 9.47 Hz, and 14.22 Hz are matched. These likely correspond to pure torsional oscillations of the HSS, which directly translate to oscillations in the generator torque. The dynamics of the IMS are not represented well with the QS method, as the QS load spectrum shows significantly lower energy in all frequencies.

In addition to the torsional oscillation peaks both the LS and KF method are able to match the gear meshing peaks. In the high-frequency range the LS method leads to a significant overestimation due to a high confidence in noisy measurements. The KF load estimates achieve a higher correlation by weighing the measurement update according to the measurement noise covariance *R* and thus filtering outliers. In the lower frequency range some peaks at 1.83 Hz, 2.91 Hz, and 6.58 Hz, which are more pronounced at the IMS-A, are missed by both the LS and KF methods.

The missed lower-frequency peaks likely relate to radial disturbance forces *f*_{dis} on the IMS, as the spectrum of measured disturbance forces suggests. The measured disturbance forces are extracted from the high-fidelity drive train simulations as connection forces of the second planetary gear stage to the IMS and show several low-frequency components of higher energy. The load estimation methods are unable to take these into account via state predictions since the disturbance forces are assumed as white gaussian process noise in the underlying physical model. Furthermore, it is challenging to consider low-frequency disturbance force excitations via vibration measurements, because these cause relatively low acceleration responses with a low signal-to-noise ratio.

For a quantitative assessment of the load correlation, the Pearson correlation coefficient is calculated for the complete time series of 3600 s, as shown in Table 2. The correlation of IMS loads is quite poor, as the studied methods are unable to reproduce aforementioned low-frequency load components. The KF is the best performing method, resulting in correlation values of 0.50–0.61. At the HSS, the QS method is sufficient to estimate bearing loads with high correlation (*R* > 0.8), as internal dynamics have less significance here. The LS and KF method do not lead to significant improvements at the HSS.

### 3.2 Fatigue Damage.

Shown in Table 3 are the relative fatigue damage errors $(D^\u2212D)/D^$ for the IMS and HSS bearings. The QS method results in low errors of 5–15% across all bearings in the considered load case. The results of the higher fidelity methods LS and KF differ only marginally from those of the QS method despite considering internal dynamics and providing load estimates of higher correlation. The error can be slightly reduced at the bearings IMS-B and HSS-B,C, however at the bearings IMS-C and HSS-A a slightly higher error is calculated. These results suggest that for the considered load case and drive train design, the fatigue damage at the IMS and HSS bearings is mainly dependent on the drive train torque and effects of internal gearbox dynamics are negligible. This becomes more clear when looking at the bearing stress cycles, which are not only a function of the load oscillations depicted in Fig. 3 but also of the rotational speed [18]. A rotating bearing experiences cyclic loading due to entering and exiting its load zone. This is reflected in the use of the stress cycle counting method LDD as opposed to the rainflow counting (RFC) method, which is commonly used for structural elements. The LDD method counts one stress cycle per revolution with a stress range equal to the current radial load. Thus, the quasi-static reactionary forces to the drive train torque cause major stress ranges and contribute significantly to the bearing fatigue, whereas the load variations from internal dynamics cause comparatively small stress ranges. In the studied load case at rated wind speed under normal operational conditions, the QS method would be sufficient to monitor fatigue damage with high accuracy and computational speed. However, it is uncertain how the QS method would perform in load cases with greater internal dynamics, such as an emergency stop or gear faults. Further studies are planned to address this topic.

## 4 Conclusion

In this article a novel approach for the estimation of wind turbine gearbox loads with the purpose of online fatigue damage monitoring was presented. The proposed method employs a Digital Twin framework and aims at continuous estimation of the dynamic states based on CMS vibration data and generator torque measurements from the SCADA data. The proposed method was evaluated in a load case at rated wind speed under normal operational conditions. With a quasi-static approach, which assumes proportionality to the drive train torque, the overall level of bearing loads was estimated with high accuracy, however the dynamic behavior was not reflected well. The quasi-static method was sufficient to estimate fatigue damage with an error of 5–15% across all bearings. The Kalman filter approach produced the highest correlation of bearing loads ranging from 0.5 to 0.96 and was able to capture high-frequency dynamics accurately, but missed several low-frequency components. These are caused by disturbance forces on the IMS, which are not reflected in the underlying physical model and are not available through measurements. Despite considering internal dynamics, the KF method did not result in significant improvements with reference to fatigue damage. It appears that in this load case the stress cycles caused by internal dynamics are insignificant relative to torque induced stress cycles. The least-squares estimator performed worse than the Kalman filter, as it is more sensitive to measurement noise. Further studies are planned to extend this work to different load cases or fault cases, assess the sensitivity to measurement noise and model uncertainties and quantify computational costs.

## Acknowledgment

The authors wish to acknowledge financial support from the Research Council of Norway through InteDiag-WTCP project (Project No. 309205).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.