On the modeling errors of digital twins for load monitoring and fatigue assessment in wind turbine drivetrains

This article presents a systematic assessment of the modeling and estimation errors of digital twins for load and fatigue monitoring in wind turbine drivetrains. The errors in the measurement input, the reduced-order drivetrain models, and the model updating methods are investigated. A statistical analysis is conducted on gear and bearing load measurements from numerical studies with 5 and 10 MW drivetrain models and from field measurements of a 1.5 MW research turbine. The error distributions are quantified using normal distributions, and limitations of the digital twin are discussed such as the information loss of 10 min averaged supervisory control and data acquisition system (SCADA) data, the estimation errors of the unknown rotor torque, and the modeling errors in torsional reduced-order drivetrain models. This study contributes to a deeper understanding of the origin and the effects of uncertainty in digital twins and delivers a foundation for further reliability and risk assessment studies.

1 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 (Wind Europe, 2020). However, major economic limitations of offshore wind turbines are high operational and maintenance expenditures (OPEX), which amount to about 34 % of the levelized cost of energy (LCOE) (Stehly and Beiter, 2020). 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 the OPEX is the geared drivetrain, with frequent failures and long downtimes, which is thus the subject of current research (, ). Digital twin (DT) is an emerging technology with prospects of decreasing the OPEX and improving the market competitiveness of offshore wind farms. The wind turbine drivetrain DT proposed by the authors in Mehlan et al. (2022a) would enable monitoring drivetrain loads and fatigue damage at otherwise inaccessible locations such as bearing and gear contacts using “virtual sensors”. A DT framework with the three components DT data, DT model, and DT decision support is envisioned for this objective (Fig. 1). The DT data comprise continuous data streams provided by the supervisory control and data acquisition system (SCADA) and the condition monitoring system (CMS), the data history including the load history and the accumulated fatigue damage, asset information such as the drivetrain topology, and general domain knowledge on drivetrain physics. The DT model refers to physics-based models to simulate internal drivetrain dynamics. Reduced-order models (ROMs) are derived from high-fidelity multibody simulation (MBS) models that are considered full-order models (FOMs) for the purpose of real-time simulation. The virtual model and its physical counterpart are synchronized with real-time field measurements using model updating techniques. State estimators such as Kalman filters are applied to infer the dynamic states of the drivetrain at small time intervals, given by the sensor sample frequency of 200 Hz. System identification methods are used to estimate system parameters such as inertia, stiffness, and damping parameters as a means to validate values provided by gearbox manufacturers or to track long-term parameter variations due to faults, material degradation, or other mechanisms. Therefore it is sufficient to update the model parameters at longer time intervals, here set to 10 min. The model updating, also referred to as data fusion or digital twinning, is essential as it facilitates the use 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 DT decision support. The focus lies on long-term fatigue damage and remaining useful life (RUL) assessment of drivetrain components, which is necessary to advance from corrective to predictive maintenance strategies. In previous numerical and field studies the proof of concept of the DT framework could be demonstrated (Mehlan et al., 2022a, 2023); however, research questions remain on the sources and the magnitude of the virtual measurements' uncertainty. Aleatory uncertainty is present in the DT data input due to the stochastic nature of wind and wave loads, and epistemic uncertainty occurs in the load and fatigue calculations due to the limitations of the DT model. Uncertainty quantification is a crucial step in the development of DTs to ensure accurate and reliable model predictions and enable informed decision-making under uncertainty (Thelen et al., 2023). The uncertainty in long-term fatigue damage calculation of wind turbine drivetrains is addressed in several studies on reliability-based design (Nejad et al., 2014; Li et al., 2017; Dong et al., 2020). Nejad et al. (2014) present a method for fatigue analysis for gear tooth root bending and differentiate between the uncertainty in the aeroelastic model, the drivetrain model, and the fatigue damage model (Nejad et al., 2014). The uncertainty is characterized by log-normal distributions with standard deviation values ranging from 0.01 for the drivetrain model to 0.1 for the aeroelastic model. Li et al. (2017) present a study on reliability-based design optimization of gear profiles and consider the uncertainty of the wind conditions with a joint probability density function of the wind speed and turbulence intensity (Li et al., 2017). Dong et al. (2020) further consider model uncertainties in a wide range of drivetrain and fatigue model parameters (Dong et al., 2020). The aforementioned studies are focused on the design of wind turbine drivetrains, where the aleatory uncertainty in the unknown environmental conditions is most influential. For DTs of operating wind turbines the challenge shifts from aleatory uncertainty towards epistemic uncertainty, since the environmental conditions and the dynamic system response are continuously estimated using real-time measurements and state estimation methods. The epistemic uncertainty of such methods has not yet been investigated systematically, as this approach is relatively novel in the field of wind energy. One investigation of the accuracy of DT-based fatigue damage monitoring is presented by Branlard et al. (2024), where errors in the range of 10 %–15 % are reported; however, the focus lies on the tower rather than drivetrain components. This study presents a numerical investigation of the epistemic uncertainty of drivetrain DTs. Simulation results from high-fidelity MBS models serve as a proxy for the real drivetrain behavior, which are developed according to best practices for component-level (gears and bearings) dynamic load calculation (Nejad et al., 2016; Wang et al., 2020). Therefore, it should be kept in mind that the reported results represent the information loss and modeling errors in relation to the high-fidelity models and only approximate the epistemic uncertainty expected in the field. Nonetheless, this analysis yields a better understanding of the sources and characteristics of epistemic uncertainty and identifies areas where optimization is possible. The remainder of this article is structured as follows: Sect. 2 presents the methodology of the DT framework for fatigue damage monitoring and defines the numerical and experimental case studies for the assessment of modeling and estimation errors. Section 3 discusses the observed error distributions in different DT components and their impact on long-term fatigue damage. Concluding remarks are provided in Sect. 4.

Figure 1 Digital twin framework for continuous remaining useful life estimation in wind turbine drivetrain components and sources of modeling and estimation errors (Mehlan et al., 2022a).

2 Methodology

2.1 Definition of modeling and estimation errors

The proposed DT framework comprises several interacting models and data processing algorithms, each of which introduces errors. These errors are grouped into the categories of measurement, state estimation, system identification, model, and fatigue damage errors. The measurement error e_meas reflects the information loss from the low sampling frequency of 10 min SCADA data, which is insufficient to observe high-frequency drivetrain dynamics. It is defined as the relative error of the measured generator torque ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{Gen}}$ to the true generator torque T_Gen. The true generator torque is sampled from simulation measurements at 200 Hz, while the measured generator torque is obtained by averaging the simulation measurements in 1 s or 10 min intervals.

$$\begin{array}{}\text{(1)}& e_{\mathrm{meas}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{T}}_{\mathrm{Gen}}-T_{\mathrm{Gen}}}{T_{\mathrm{Gen}}}}\end{array}$$

The state estimation error e_SE refers to errors caused by the Kalman filter algorithm. The Kalman filter fuses uncertain information from measurements and model predictions and is the optimal state estimator in the case of white Gaussian measurement and process noise. However, this assumption is not valid here since the unknown rotor torque modeled as process noise exhibits nonuniformity such as peaks at characteristic excitation frequencies (1P, 3P, …). It is therefore expected that use of a Kalman filter introduces additional errors in the drivetrain case. This error is defined as the relative error of the estimated states $\widehat{x}$ to the true states x.

$$\begin{array}{}\text{(2)}& e_{\mathrm{SE}}={\displaystyle \frac{\widehat{x}-x}{x}}\end{array}$$

The system identification error e_SI reflects the error that is introduced by the inverse methods to estimate the system's inertia, stiffness, and damping matrices ($\widehat{J},\widehat{K},\widehat{C}$) and is defined as the relative error of the estimated parameters $\widehat{{\mathit{\theta }}_i}$ to the true parameters θ_i.

$$\begin{array}{}\text{(3)}& e_{\mathrm{SI},i}={\displaystyle \frac{{\widehat{\mathit{\theta }}}_i-{\mathit{\theta }}_i}{{\mathit{\theta }}_i}}\end{array}$$

The modeling errors e_model,j refer to the limitations of the DT model to simulate all relevant drivetrain dynamics. ROMs with a limited number of torsional degrees of freedom (DOFs) are considered, which are unable to capture non-torsional drivetrain dynamics such as shaft bending modes or complex torsional dynamics such as gear meshing. The error caused by the model complexity reduction is described with the modeling errors and defined as the relative error of the ROM loads ${\widehat{F}}_j$ to the FOM loads F_j in the drivetrain components j.

$$\begin{array}{}\text{(4)}& e_{\mathrm{model},j}={\displaystyle \frac{{\widehat{F}}_j-F_j}{F_j}}\end{array}$$

Lastly, the propagation of errors through the DT's computational chain are quantified with the relative error in long-term fatigue damage e_D,j for each drivetrain component j.

$$\begin{array}{}\text{(5)}& e_{\mathrm{D},j}={\displaystyle \frac{{\widehat{D}}_j-D_j}{D_j}}\end{array}$$

The estimated fatigue damage ${\widehat{D}}_j$ includes errors from the measurement input, the state estimation method, the system identification method, and the ROM modeling errors, while the true fatigue damage is based on FOM simulations. Errors in the fatigue damage model itself such as errors introduced by the stress cycle counting method and the S–N curves are not considered here.

2.2 Numerical case studies

Two numerical case studies with the National Renewable Energy Laboratory (NREL) 5 MW baseline turbine (Jonkman et al., 2009) and the DTU 10 MW reference turbine are conducted (Bak et al., 2013). The best practice for dynamic drivetrain simulation is the decoupled analysis approach, where the “global” structural blade and tower dynamics and the “internal” drivetrain dynamics are simulated separately (Nejad et al., 2014). The global system response is simulated first with an aeroelastic model, and the resultant main shaft loads and nacelle motions are then imposed as boundary conditions on the drivetrain model. This procedure is motivated by the fact that the global dynamics are governed by aerodynamic excitations and occur at low frequencies (<10 Hz), while much higher frequencies such as gear meshing frequencies at >100 Hz need to be considered for the drivetrain dynamics. The simulation cases are designed according to the IEC 61400-1 requirements for long-term fatigue analysis. A total of 12 wind speed cases ranging from cut-in wind speed of 3 m s⁻¹ to cut-out wind speed of 25 m s⁻¹ are considered (Table 1). One case of turbulence intensity is considered and modeled with the IEC turbulence class A. Only one case of wave height and wave period is considered, since the drivetrain bearing and gear loads are reportedly insensitive to the sea state. The primary effect of harsher sea states can be observed in increased axial loads induced by pitch motions, which are compensated for by the main bearings in a four-point suspension and do not propagate further into the drivetrain (Nejad et al., 2015). Each environmental condition (EC) is simulated for 1 h with six different random realizations (seeds) of turbulent wind fields.

Table 1 Environmental conditions for simulation with global and drivetrain models.

2.3 Global models

The global wind turbine dynamics are simulated with open-source aeroelastic models of the NREL 5 MW and DTU 10 MW reference turbines mounted on the OC4 and Nautilus semisubmersible platforms, respectively (Robertson et al., 2012; Arias and Galvan, 2018). The models are implemented in the aeroelastic code OpenFAST that comprises computational modules for calculation of the aerodynamics, hydrodynamics, structural dynamics, and wind turbine control (OpenFAST, 2022). The aerodynamics are calculated with blade element momentum (BEM) theory, where the turbulent wind field is generated with the Kaimal turbulence model. The structural dynamics of the blades and the tower are based on Timoshenko elastic beam theory. The incident wave loads on the floater are modeled with a Jonswap spectrum. A variable-speed controller is implemented for the 5 MW and 10 MW model.

2.4 Full order drivetrain models

Multibody simulation (MBS) models of the NREL 5 MW and DTU 10 MW reference turbine serve as a benchmark in this study (Nejad et al., 2016; Wang et al., 2020). The MBS models are developed according to best practices and current model fidelity guidelines (Guo et al., 2015) and are thus considered FOMs. Both FOMs have similar topology and comprise a four-point suspension for the main shaft and a gearbox with two planetary gear stages and one parallel gear stage (Fig. 2). However, the 5 MW model represents a high-speed gearbox with a gear ratio of 1:96.354, while the 10 MW model represents a medium-speed gearbox with a gear ratio of 1:50.039. The FOMs allow shaft motion in all six degrees of freedom (DOFs) and consider the flexibility in the main shaft and the planet carriers. The flexible bodies are implemented as condensed FE models with modal damping of 2 % (Guo et al., 2015). The bearings and the torque arm bushings are modeled as linear spring–damper connections in six DOFs with diagonal stiffness and damping matrices, where the damping is assumed to be stiffness-proportional in the range of 0.25 % to 2.5 % (Nejad et al., 2016). The gear compliance is modeled with a time-invariant mesh stiffness function capable of emulating gear meshing excitations. Stiffness-proportional damping with a factor of 0.1 % is added (Mehlan et al., 2022b). The input loads simulated with aeroelastic models are imposed on the main shaft, while the generator shaft speed is controlled with a PI controller.

Figure 2 Topology and component nomenclature of the NREL 5 MW and DTU 10 MW drivetrain models.

2.5 Reduced-order drivetrain models

Reduced-order models (ROMs) are preferable as DT models due to the high computational costs in real-time monitoring applications (Mehlan et al., 2022b). The complexity of DT models is also limited by the observability requirement of the state estimator. The state estimator that is used to match the dynamics of the DT model with the physical turbine requires that all dynamic states are observable with the available measurement input. The SCADA measurements of the main and generator shaft speeds allow the observation of torsional drivetrain modes. Bending and lateral drivetrain modes are observable with CMS accelerometers mounted on the gearbox housing; however, the sensitivity is relatively low due to measurement noise and the observation function is complex due to the transfer path of the vibration through the housing (Mehlan et al., 2022b). For this reason, the ROMs are limited to torsional degrees of freedom (DOFs) only. Lumped parameter models with one and two torsional DOFs are considered. The input torques at each gear stage T_in,k are calculated with the torsional ROMs and then further used to determine local gear and bearing forces (Sect. 2.5.3).

2.5.1 Rigid one-degree-of-freedom ROM

The first ROM represents a rigid, torsional model with one degree of freedom (DOF). The flexibility of shafts and gear contacts are neglected, which yields direct coupling of the angular shaft velocities ω_k and input torques at each gear stage T_in,k via the gear ratios i_k.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{\omega }}_{\mathrm{Rot}}& ={\mathit{\omega }}_{\mathrm{in},\mathrm{2}}/i_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{in},\mathrm{3}}/i_{\mathrm{1}}/i_{\mathrm{2}}={\mathit{\omega }}_{\mathrm{Gen}}/i_{\mathrm{1}}/i_{\mathrm{2}}/i_{\mathrm{3}}\\ T_{\mathrm{in},\mathrm{1}}& =i_{\mathrm{1}}{T}_{\mathrm{in},\mathrm{2}}=i_{\mathrm{1}}{i}_{\mathrm{2}}{T}_{\mathrm{in},\mathrm{3}}=i_{\mathrm{1}}{i}_{\mathrm{2}}{i}_{\mathrm{3}}{T}_{\mathrm{Gen}}\end{array}\end{array}$$

The rigid ROM is advantageous, in that it does not require inertia, stiffness, or damping parameters for model construction and validation, which minimizes the errors associated with system identification techniques for parameter estimation (e_SI). In addition, it is not necessary to apply state estimation methods, since the gear stage torques and thus all drivetrain loads are directly observable with the measured generator torque, which reduces state estimation errors (e_SE).

2.5.2 Flexible two-degrees-of-freedom ROM

The second ROM introduces one additional torsional DOF and is able to represent the first torsional mode. However, this model assumes knowledge of inertia, stiffness, and damping parameters, which may be estimated via system identification techniques. The flexibility of all drivetrain components is lumped into a scalar drivetrain stiffness k_DT, while the torsional inertias are lumped into either the rotor inertia J_Rot or the generator inertia J_Gen. The equations of motion are then given by

$$\begin{array}{}\text{(7)}& \mathbf{J}\ddot{\mathit{\varphi }}+\mathbf{C}\dot{\mathit{\varphi }}+\mathbf{K}\mathit{\varphi }+\mathit{f}=\mathrm{0},\end{array}$$

where J denotes the inertia matrix, C is the damping matrix, K is the stiffness matrix, f is the external force vector, and ϕ represents the independent dynamic states:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathbf{J}=\left[\begin{array}{cc}{J}_{\mathrm{Rot}}& \mathrm{0}\\ \mathrm{0}& J_{\mathrm{Gen}}\end{array}\right],\mathbf{C}=\left[\begin{array}{cc}{c}_{\mathrm{DT}}& -c_{\mathrm{DT}}/i_{\mathrm{DT}}\\ -c_{\mathrm{DT}}/i_{\mathrm{DT}}& c_{\mathrm{DT}}/i_{\mathrm{DT}}^{\mathrm{2}}\end{array}\right],\\ & \mathbf{K}=\left[\begin{array}{cc}{k}_{\mathrm{DT}}& -k_{\mathrm{DT}}/i_{\mathrm{DT}}\\ -k_{\mathrm{DT}}/i_{\mathrm{DT}}& k_{\mathrm{DT}}/i_{\mathrm{DT}}^{\mathrm{2}}\end{array}\right],\\ & \mathit{\varphi }=\left[\begin{array}{c}{\mathit{\varphi }}_{\mathrm{Rot}}\\ {\mathit{\varphi }}_{\mathrm{Gen}}\end{array}\right],\mathit{f}=\left[\begin{array}{c}-T_{\mathrm{Rot}}\\ T_{\mathrm{Gen}}\end{array}\right].\end{array}\end{array}$$

The gear stage input torques are still coupled and only a function of the rotor and generator shaft angular positions ϕ.

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{T}_{\mathrm{in},\mathrm{1}}& =i_{\mathrm{1}}{T}_{\mathrm{in},\mathrm{2}}=i_{\mathrm{1}}{i}_{\mathrm{2}}{T}_{\mathrm{in},\mathrm{3}}=[c_{\mathrm{DT}},-c_{\mathrm{DT}}/i_{\mathrm{DT}}]\dot{\mathit{\varphi }}\\ & +[k_{\mathrm{DT}},-k_{\mathrm{DT}}/i_{\mathrm{DT}}]\mathit{\varphi }\end{array}\end{array}$$

2.5.3 Bearing and gear forces

The gear forces are determined with free-body diagrams and moment balances as a function of the gear stage input torques. Dynamic effects of planet load sharing are not considered at the planetary gear stages; hence, the gear stage torque is distributed equally among the number of planets N_PL. Furthermore, the gear forces at the ring–planet and the sun–planet contacts are assumed to be equal. The circumferential (z-direction) gear forces F_t are then obtained as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{F}_{t,\mathrm{1}}& =T_{\mathrm{in},\mathrm{1}}\cdot i_{\mathrm{1}}/r_{\mathrm{b},\mathrm{S},\mathrm{1}}/N_{\mathrm{PL},\mathrm{1}}\\ F_{t,\mathrm{2}}& =T_{\mathrm{in},\mathrm{2}}\cdot i_{\mathrm{2}}/r_{\mathrm{b},\mathrm{S},\mathrm{2}}/N_{\mathrm{PL},\mathrm{2}}\\ F_{t,\mathrm{3}}& =T_{\mathrm{in},\mathrm{3}}/r_{\mathrm{b},\mathrm{G},\mathrm{3}},\end{array}\end{array}$$

where r_b represents the base radii of the first- and second-stage sun and of the third-stage gear wheel. The remaining gear force components in the x and y direction, the axial and radial gear force components F_a and F_r, are determined with the tangential pressure angle α_t and helix angle β. The planetary gear stage is modeled with spur gears (β=0), while the parallel gear stage is modeled with a helix angle of β=10°.

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{F}_{\mathrm{r}}& =F_t\tan \left({\mathit{\alpha }}_{\mathrm{t}}\right)/\cos \left(\mathit{\beta }\right)\\ F_{\mathrm{a}}& =F_t\tan \left(\mathit{\beta }\right)\end{array}\end{array}$$

At the planetary gear stages the radial bearing forces F_rad are directly proportional to the circumferential gear forces with the assumption of negligible gravity forces.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& F_{\text{rad, PL-A}}=\mathrm{2}\cdot F_{t,\mathrm{1}}\\ & F_{\text{rad, IMS-PL-A}}=\mathrm{2}\cdot F_{t,\mathrm{2}}\end{array}\end{array}$$

At the helical gear stage the radial bearing forces are derived with moment balances:

$$\begin{array}{}\text{(13)}& F_{\mathrm{rad}}=\sqrt{F_y^{\mathrm{2}}+F_z^{\mathrm{2}}},\end{array}$$

with the following.

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}& F_{y,\text{IMS-A}}=-F_{\mathrm{r}}{\displaystyle \frac{d_{\text{IMS-B}}-d_{\text{W}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}+F_{\mathrm{a}}{\displaystyle \frac{r_{\text{p,W}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}\\ & F_{y,\text{IMS-B}}=-F_{\mathrm{r}}{\displaystyle \frac{d_{\text{W}}-d_{\text{IMS-A}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}-F_{\mathrm{a}}{\displaystyle \frac{r_{\text{p,W}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}\\ & F_{y,\text{HSS-A}}=F_{\mathrm{r}}{\displaystyle \frac{d_{\text{HSS-B}}-d_{\text{P}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}+F_{\mathrm{a}}{\displaystyle \frac{r_{\text{p,P}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}\\ & F_{y,\text{HSS-B}}=F_{\mathrm{r}}{\displaystyle \frac{d_{\text{P}}-d_{\text{HSS-A}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}-F_{\mathrm{a}}{\displaystyle \frac{r_{\text{p,P}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}\end{array}\text{(15)}& {\displaystyle }\begin{array}{rl}& F_{z,\text{IMS-A}}=-F_{t,\mathrm{3}}{\displaystyle \frac{d_{\text{IMS-B}}-d_{\text{W}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}\\ & F_{z,\text{IMS-B}}=-F_{t,\mathrm{3}}{\displaystyle \frac{d_{\text{W}}-d_{\text{IMS-A}}}{d_{\text{IMS-B}}-d_{\text{IMS-A}}}}\\ & F_{z,\text{HSS-A}}=F_{t,\mathrm{3}}{\displaystyle \frac{d_{\text{HSS-B}}-d_{\text{P}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}\\ & F_{z,\text{HSS-B}}=F_{t,\mathrm{3}}{\displaystyle \frac{d_{\text{P}}-d_{\text{HSS-A}}}{d_{\text{HSS-B}}-d_{\text{HSS-A}}}}\end{array}\end{array}$$

The axial gear force component of the helical high-speed gear stage is supported by the HSS-B and IMS-B bearings.

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& F_{\text{ax, IMS-B}}=F_{\mathrm{a}}\\ & F_{\text{ax, HSS-B}}=-F_{\mathrm{a}}\end{array}\end{array}$$

2.6 Experimental case study

The simulation measurements are partially validated with field measurements of the Department of Energy (DOE) 1.5 MW research turbine located at the National Renewable Energy Laboratory (NREL) (Santos and van Dam, 2015). The DOE 1.5 MW turbine is equipped with a commercial Winergy PEAB 4410.4 high-speed gearbox with similar three-stage topology as the above simulation models. The dataset, originally collected for the analysis of cage and roller slip in the HSS-A bearing (Guo and Keller, 2020), is repurposed for this study. The original sample frequency of 5 kHz necessary to observe slip dynamics restricted the measurement duration and as a result the total recorded data amount to only about 30 min. Nonetheless, the amount of data for each wind speed bin is sufficient to characterize the dynamic behavior in the full range of operational conditions, which allows the comparison with simulated data from the numerical case studies. The loading of the HSS is fully determined with three shaft-mounted strain gauge bridges, one for measuring torque and two 90° offset bridges for measuring bending moments. The forces at the HSS-A bearing are calculated with the torque T and bending moment measurements M_y and M_z (Guo and Keller, 2020).

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& F_{y,\text{HSS-A}}={\displaystyle \frac{\mathrm{1}}{d_{\text{B}}-d_{\text{HSS-A}}}}\left[-M_z-T/r_{\mathrm{b}}(d_{\text{B}}-d_{\text{P}})\sin \mathit{\beta }\right]\\ & F_{z,\text{HSS-A}}={\displaystyle \frac{\mathrm{1}}{d_{\text{B}}-d_{\text{HSS-A}}}}\left[-M_y-T/r_{\mathrm{b}}(d_{\text{B}}-d_{\text{P}})\cos \mathit{\beta }\right]\end{array}\end{array}$$

These measurements are considered FOM bearing load measurements, since all relevant torsional and shaft bending dynamics are captured. The FOM loads are set in relation to the ROM loads calculated solely with torque measurements and the rigid ROM (Sect. 2.5.1) to assess the modeling errors.

Figure 3 Forces at the HSS.

2.7 State and input estimation

The DT model is synchronized with the operating wind turbine at regular time intervals Δt such that gear and bearing loads can be measured with “virtual sensors” in the synchronized model. The challenge lies in the incomplete and noisy measurements of both the dynamic states and the input forces, which poses a joint state and input estimation problem. The measurements of the dynamic states, the shaft angular velocities and positions, are corrupted with measurement noise, while the input forces at the main shaft are unknown; only the generator side torque is measured. The augmented Kalman filter is applied here as a joint state and input estimator, as it is the optimal estimator for dynamic systems governed by linear, stochastic equations subjected to white Gaussian process and measurement noise. For this purpose the equations of motion of the flexible ROMs are first brought into discrete state-space representation:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathit{x}}_{n+\mathrm{1}}={\mathbf{F}}^{\mathrm{d}}{\mathit{x}}_n+{\mathbf{G}}_{\mathrm{k}}^{\mathrm{d}}{\mathit{u}}_{k,n}+{\mathbf{G}}_{\mathrm{u}}^{\mathrm{d}}{\mathit{u}}_{\mathrm{u},n}+{\mathit{w}}_n,\text{(19)}& {\displaystyle }{\mathit{y}}_n={\mathbf{H}}^{\mathrm{d}}{\mathit{x}}_n+{\mathit{v}}_n,\end{array}$$

where the state vector x is obtained by stacking the shaft angular positions and velocities, the input forces u are split into the known generator torque u_k and the unknown rotor torque u_u, the measurement vector y contains the rotor and generator shaft speeds, the unknown dynamic component of the rotor torque is considered white Gaussian process noise w with covariance Q, and v is white Gaussian measurement noise with covariance R.

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& \mathit{x}=[\mathit{\varphi }\dot{\mathit{\varphi }}{]}^{\mathrm{T}}\\ & {\mathit{u}}_{\mathrm{k}}=T_{\mathrm{Gen}}\\ & {\mathit{u}}_{\mathrm{u}}=T_{\mathrm{Rot}}\\ & \mathit{y}=[{\dot{\mathit{\varphi }}}_{\mathrm{Rot}},{\dot{\mathit{\varphi }}}_{\mathrm{Gen}}{]}^{\mathrm{T}}\\ & \mathit{w}\sim \mathcal{N}(\mathbf{0},\mathbf{Q})\\ & \mathit{v}\sim \mathcal{N}(\mathbf{0},\mathbf{R})\end{array}\end{array}$$

The system matrix F^d, the input matrix G^d, and the observation matrix H^d of the discrete state-space model are calculated as follows

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathbf{F}}^{\mathrm{d}}& {\displaystyle }=\exp \left({\mathbf{F}}^{\mathrm{c}}\mathrm{\Delta }t\right),\text{(22)}& {\displaystyle }{\mathbf{G}}^{\mathrm{d}}& {\displaystyle }=\left[{\mathbf{G}}_{\mathrm{k}}^{\mathrm{d}}{\mathbf{G}}_{\mathrm{u}}^{\mathrm{d}}\right]=\left({\mathbf{F}}^{\mathrm{c}}{)}^{-\mathrm{1}}\right({\mathbf{F}}^{\mathrm{d}}-{\mathbf{I}}^{\mathrm{2}N\times \mathrm{2}N})\left[{\mathbf{G}}_{\mathrm{k}}^{\mathrm{c}}{\mathbf{G}}_{\mathrm{u}}^{\mathrm{c}}\right],\text{(23)}& {\displaystyle }{\mathbf{H}}^{\mathrm{d}}& {\displaystyle }=\left[\begin{array}{cc}{\mathbf{0}}^{N\times N}& {\mathbf{I}}^{N\times N}\end{array}\right],\end{array}$$

where N denotes the model's DOF, 0 is the null matrix, I is the identity matrix, and F^c, ${\mathbf{G}}_{\mathrm{k}}^{\mathrm{c}}$, ${\mathbf{G}}_{\mathrm{u}}^{\mathrm{c}}$, and H^c are the matrices of the continuous state-space model.

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& {\mathbf{F}}^{\mathrm{c}}=\left[\begin{array}{cc}{\mathbf{0}}^{N\times N}& {\mathbf{I}}^{N\times N}\\ -{\mathbf{J}}^{-\mathrm{1}}\mathbf{K}& -{\mathbf{J}}^{-\mathrm{1}}\mathbf{C}\end{array}\right]\\ & {\mathbf{G}}_{\mathrm{k}}^{\mathrm{c}}={\left[\begin{array}{cc}{\mathbf{0}}^{\mathrm{1}\times \mathrm{2}N-\mathrm{1}}& \mathrm{1}/J_{\mathrm{Gen}}\end{array}\right]}^{\mathrm{T}}\\ & {\mathbf{G}}_{\mathrm{u}}^{\mathrm{c}}={\left[\begin{array}{ccc}{\mathbf{0}}^{\mathrm{1}\times N}& -i_{\mathrm{DT}}/J_{\mathrm{Rot}}& {\mathrm{0}}^{\mathrm{1}\times \mathrm{N}-\mathrm{1}}\end{array}\right]}^{\mathrm{T}}\\ & {\mathbf{H}}^{\mathrm{c}}={\mathbf{H}}^{\mathrm{d}}\end{array}\end{array}$$

For the purpose of simultaneous state and input estimation, the state vector x is expanded with the unknown input force u_u, yielding the state-space representation with the augmented state vector x^a=[xu_u]^T.

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathit{x}}_{n+\mathrm{1}}^{\mathrm{a}}=\mathbf{F}{\mathit{x}}_n^{\mathrm{a}}+{\mathbf{G}}_{\mathrm{k}}{\mathit{u}}_{k,n}+{\mathit{w}}_n,\text{(26)}& {\displaystyle }{\mathit{y}}_n=\mathbf{H}{\mathit{x}}_n^{\mathrm{a}}+{\mathit{v}}_n,\end{array}$$

where the system matrix F, the input matrix G, and the observation matrix H of the augmented state-space model are calculated as follows.

$$\begin{array}{}\text{(27)}& {\displaystyle }\mathbf{F}=\left[\begin{array}{c}{\mathbf{F}}^{\mathbf{d}}{\mathbf{G}}_{\mathbf{u}}^{\mathbf{d}}\\ {\mathbf{0}}^{\mathrm{1}\times \mathrm{N}}\mathrm{1}\end{array}\right]\text{(28)}& {\displaystyle }\mathbf{G}=\left[\begin{array}{c}{\mathbf{G}}_{\mathrm{k}}^{\mathrm{d}}\\ \mathrm{0}\end{array}\right]\text{(29)}& {\displaystyle }\mathbf{H}=\left[\begin{array}{c}{\mathbf{H}}^{\mathrm{d}}{\mathbf{0}}^{N\times \mathrm{1}}\end{array}\right]\end{array}$$

The Kalman filter produces the state estimates $\widehat{\mathit{x}}$ in a two-step algorithm, comprising the prediction step and the measurement update step.

$$\begin{array}{}\text{(30)}& {\displaystyle }{\widehat{\mathit{x}}}_{n|n-\mathrm{1}}^{\mathrm{a}}=\mathbf{F}{\widehat{\mathit{x}}}_{n-\mathrm{1}|n-\mathrm{1}}^{\mathrm{a}}+\mathbf{G}{\mathit{u}}_{n-\mathrm{1}},\text{(31)}& {\displaystyle }{\widehat{\mathbf{P}}}_{n|n-\mathrm{1}}=\mathbf{F}{\widehat{\mathit{x}}}_{n-\mathrm{1}|n-\mathrm{1}}^{\mathrm{a}}{\mathbf{F}}^{\mathrm{T}}+\mathbf{Q},\text{(32)}& {\displaystyle }{\mathbf{M}}_n={\widehat{\mathbf{P}}}_{n|n-\mathrm{1}}{\mathbf{H}}^{\mathrm{T}}(\mathbf{H}{\widehat{\mathbf{P}}}_{n|n-\mathrm{1}}{\mathbf{H}}^{\mathrm{T}}+\mathbf{R}{)}^{-\mathrm{1}},\text{(33)}& {\displaystyle }{\widehat{\mathit{x}}}_{n|n}^{\mathrm{a}}={\widehat{\mathit{x}}}_{n|n-\mathrm{1}}^{\mathrm{a}}+{\mathbf{M}}_n({\mathit{y}}_n-\mathbf{H}{\widehat{\mathit{x}}}_{n|n-\mathrm{1}}^{\mathrm{a}}),\text{(34)}& {\displaystyle }{\widehat{\mathbf{P}}}_{\mathbf{n}|\mathbf{n}}=(\mathbf{I}-{\mathbf{M}}_n\mathbf{H}){\widehat{\mathbf{P}}}_{n|n-\mathrm{1}}.\end{array}$$

2.8 System identification

System identification methods are applied to continuously update the model properties to ensure the convergence of the virtual model and the physical wind turbine's dynamic behavior. The rotor inertia, generator inertia, drivetrain torsional stiffness, and damping are considered time-variant parameters to reflect long-term changes in the physical wind turbine. The rotor inertia may increase due to the accretion of dirt, moisture, and ice or decrease as a result of leading-edge erosion or similar damages. The drivetrain stiffness and damping values are affected by material fatigue and localized faults such as spalls or tooth root cracks. The second line of the equations of motion (Eq. 7) is used to estimate the parameter set $\mathit{\theta }=[J_{\mathrm{Gen}},c_{\mathrm{DT}},k_{\mathrm{DT}},{\mathit{\alpha }}_{\mathrm{0}}]$, since the boundary conditions are fully determined here by measurements of the generator torque. The following least-squares optimization problem is then formulated.

$$\begin{array}{}\text{(35)}& \begin{array}{rl}\widehat{\mathit{\theta }}& =\arg \underset{\mathit{\theta }}{min}\left|\right|J_{\mathrm{Gen}}{\ddot{\mathit{\varphi }}}_{\mathrm{Gen}}-c_{\mathrm{DT}}/i_{\mathrm{DT}}\dot{\mathit{\alpha }}\\ & -k_{\mathrm{DT}}/i_{\mathrm{DT}}(\mathit{\alpha }-{\mathit{\alpha }}_{\mathrm{0}})+T_{\mathrm{Gen}}|{|}_{\mathrm{2}}^{\mathrm{2}}\end{array}\end{array}$$

The generator shaft acceleration ${\ddot{\mathit{\varphi }}}_{\mathrm{Gen}}$ is obtained by numerical differentiation of the measured SCADA generator shaft speed. The drivetrain torsion defined as $\mathit{\alpha }={\mathit{\varphi }}_{\mathrm{Rot}}-{\mathit{\varphi }}_{\mathrm{Gen}}/i_{\mathrm{DT}}$ is calculated by numerical integration of the shaft speeds. As a result of the numerical integration of noisy signals, a runaway trend or sensor drift is observed, which is removed via MATLAB's detrend function. Furthermore, the initial state α₀ of the integrated signal is unknown and therefore added to the parameter set of the optimization problem. The optimization problem is solved for 10 min time sections at each EC using a least-squares solver. Unfortunately, the same procedure cannot be employed to obtain the remaining parameter, the rotor inertia J_Rot, since the rotor torque is typically not measured by SCADA systems, which leaves the rotor side equations of motion undefined (Eq. 7). Operational modal analysis (OMA) techniques are used instead. The first torsional natural frequency ${\widehat{f}}_N$ is estimated using peak finding algorithms in the frequency spectrum of the drivetrain torsion signal α. Since the natural frequency is a function of the drivetrain inertia and stiffness, one may solve for the unknown rotor inertia as follows.

$$\begin{array}{}\text{(36)}& \begin{array}{rl}& {\widehat{J}}_{\mathrm{eq}}={\displaystyle \frac{{\widehat{k}}_{\mathrm{DT}}}{(\mathrm{2}\mathit{\pi }{\widehat{f}}_N{)}^{\mathrm{2}}}}\\ & {\widehat{J}}_{\mathrm{Rot}}=(\mathrm{1}/{\widehat{J}}_{\mathrm{eq}}-\mathrm{1}/{\widehat{J}}_{\mathrm{Gen}}/i_{\mathrm{DT}}{)}^{-\mathrm{1}}\end{array}\end{array}$$

2.9 Fatigue damage

The gear and bearing fatigue damage is based on the gear tooth root stress calculation of ISO 6336 (ISO 6336, 2006) and the nominal bearing life calculation of ISO 281 (ISO 281, 2007). The gear tooth root stress s is determined from the circumferential gear force F_t, the flank width b, the normal module m_n, and the modification factors Y and K (ISO 6336, 2006).

$$\begin{array}{}\text{(37)}& s={\displaystyle \frac{F_{\mathrm{t}}}{b{m}_{\mathrm{n}}}}{Y}_{\mathrm{S}}{Y}_{\mathrm{F}}{Y}_{\mathit{\beta }}{Y}_{\mathrm{B}}{Y}_{\mathrm{DT}}{K}_{\mathrm{A}}{K}_{\mathrm{V}}{K}_{\mathrm{F}\mathit{\beta }}{K}_{{\mathrm{F}}_{\mathit{\alpha }}}{K}_{\mathit{\gamma }}\end{array}$$

The pendant for bearings is the equivalent dynamic load P that is defined for cylindrical roller bearings (CRBs) and tapered roller bearings (TRBs) as follows (ISO 281, 2007).

$$\begin{array}{}\text{(38)}& {\displaystyle }\text{CRB:}P=F_{\mathrm{rad}},\text{(39)}& {\displaystyle }\text{TRB:}P=\left\{\begin{array}{l}{F}_{\mathrm{rad}}+Y_{\mathrm{1}}{F}_{\mathrm{ax}},\text{if }{F}_{\mathrm{ax}}/F_{\mathrm{rad}}\le e\\ \mathrm{0.67}{F}_{\mathrm{rad}}+Y_{\mathrm{2}}{F}_{\mathrm{ax}},\text{otherwise},\end{array}\right.\end{array}$$

where Y₁, Y₂, and e are bearing-specific parameters. The load duration distribution (LDD) method is used as a stress cycle counting method for components in rotating machinery that experience cyclic loading due to entering and exiting the load zone (Nejad et al., 2014). The LDD method counts one stress cycle per shaft revolution and distributes the cycles n_i into 64 bins of increasing stress range. The permissible stress cycles N_i for each stress range are modeled with S–N curves for gear tooth root fatigue,

$$\begin{array}{}\text{(40)}& N_i=K_{\mathrm{c}}{s}_i^{-m},\end{array}$$

where m=6.225 and K_c=10^24.744 (Nejad et al., 2014), and the nominal bearing life equation for bearing fatigue (ISO 281, 2007),

$$\begin{array}{}\text{(41)}& N_i={\mathrm{10}}^{\mathrm{6}}{\left({\displaystyle \frac{C}{P_i}}\right)}^m,\end{array}$$

where C is the basic dynamic load rating and $m=\mathrm{10}/\mathrm{3}$ for roller bearings. The short-term fatigue damage is then calculated for 10 min time sections by summation of all stress range bins.

$$\begin{array}{}\text{(42)}& D^{\mathrm{ST}}=\sum _i{n}_i/N_i\end{array}$$

The long-term fatigue damage D^LT for the nominal lifetime of 20 years is extrapolated from the short-term fatigue damage by weighting with the wind speed probability density function f(u_k). A representative wind speed distribution measured at Anholt, Denmark, is selected (Gonzalez et al., 2019b).

$$\begin{array}{}\text{(43)}& D^{\mathrm{LT}}={\displaystyle \frac{\mathrm{20}\mathrm{year}}{\mathrm{10}\min }}\sum _{\mathrm{k}}f\left(u_{\mathrm{k}}\right)D_{\mathrm{k}}^{\mathrm{ST}}\end{array}$$

3 Results and discussions

3.1 Choice of error distribution

The first step in the statistical analysis is the identification of the error distribution shapes, which are of importance in reliability and risk assessment studies. A common assumption for dynamic modeling errors and estimation errors in Kalman filters is to use Gaussian distributions according to the central limit theorem. To check this assumption, the Gaussian distribution is benchmarked against 13 different commonly used statistical distributions. The fitted distributions for the measurement, state estimation, system identification, and modeling errors are ranked according to their goodness of fit given by the coefficient of determination R². Figure 4 shows the R² values of the six best-performing distributions aggregated for all ECs of the 5 MW case study. The results are inconclusive as to which distribution is best suited, but it can be stated that the normal distribution yields a reasonable fit of R²>0.9 for all types of modeling and estimation errors in DTs. The further statistical analysis is continued with normal distributions for interpretability and comparability with other publications. The findings of this study are summarized in Table 3 and discussed in the following sections.

Figure 4 Goodness of fit of different distribution shapes for the measurement, state estimation, system identification, and modeling errors, aggregated for all ECs of the 5 MW case study.

3.2 Measurement errors

The first source of errors in the proposed load and fatigue monitoring approach originates from the low temporal resolution of the SCADA data input. Typical SCADA systems operate with sampling frequencies of 1 Hz but store the data only as 10 min averages, which has already been identified as a limiting factor for monitoring approaches. The generator torque reportedly has the fastest decaying autocorrelation out of all SCADA signals, which results in a large loss of information when using time-averaged signals (Gonzalez et al., 2019a). This motivated efforts in the industry to adopt high-frequency (1 Hz) SCADA systems; however, even a sampling frequency of 1 Hz is arguably insufficient to fully capture drivetrain dynamics, since the first torsional natural frequency and internal excitation frequencies such as gear meshing frequencies lie well above the Nyquist frequency of 0.5 Hz. The effects of this are illustrated in Fig. 5, which shows the standard deviation of the measurement error distributions e_meas resulting from either 1 s or 10 min averaging of the generator torque input. The information loss of 10 min data is particularly high below rated wind speed and reaches values of up to 74 % near cut-in wind speed. In wind turbines with variable-speed controllers this operational regime is characterized by a high variance in the drivetrain torque, which is not reflected in 10 min averaged data. The information loss of 1 Hz data amounts to less than 5 % for all operational conditions, which suggests that this resolution is sufficient to observe low-frequency (<0.5 Hz) load variations due to the wind speed volatility. The remaining errors are related to neglecting higher-frequency dynamics such as torsional drivetrain modes. Based on these results a measurement resolution of at least 1 Hz is recommended for load and fatigue damage monitoring in wind turbine drivetrains.

Figure 5 Standard deviation of measurement error distributions e_meas indicating the information loss of a 1 Hz or 10 min sampling frequency.

3.3 State estimation errors

The second source of errors is also related to the limitations of the SCADA measurements, in that the rotor torque is typically not measured and must be estimated indirectly using the augmented Kalman filter as a joint input–state estimation method. The state estimation errors are illustrated in Fig. 6, which shows the fitted standard deviation for the numerical case studies with the 5 and 10 MW model at different ECs. Higher errors of 3.5 % are observed at the cut-in wind speed, which can be attributed to start-up and shut-down effects. Above rated wind speed the errors exhibit a progressive trend and reach maximum values of 7 % at the cut-out wind speed. A slightly higher error is observed with the 5 MW model, which is also apparent in the frequency spectra and time series shown in Fig. 7. The rotor torque estimates for the 10 MW turbine show good agreement in the low-frequency range and at the peaks of the first torsional natural frequency (2.08 Hz). For the 5 MW turbine, on the other hand, the rotor torque is underestimated at the first torsional natural frequency (1.7 Hz) and at higher-order modes.

Figure 6 Standard deviation of state estimation error distributions e_SE.

Figure 7 True and estimated rotor torque (T_rot and ${\widehat{T}}_{\mathrm{rot}}$) using joint state–input estimation methods. Shown are the PSD frequency spectrum and the time series at EC8.

3.4 System identification errors

The third source of errors originates from the inaccuracies in the system identification method that is employed to update the DT's model parameters ($J_{\mathrm{Rot}},J_{\mathrm{Gen}},k_{\mathrm{DT}},c_{\mathrm{DT}}$). The system identification error is quantified by comparison of the parameter estimates with reference values of the 5 and 10 MW FOMs. The FOMs comprise many components such as flexible shafts, gears, and bearings that contribute in different degrees to the total inertia, stiffness, and damping (Sect. 2.4). Therefore, to obtain comparable reference values, the equivalent rotor inertia, generator inertia, drivetrain torsional stiffness, and damping are adopted from Wang et al. (2020) and Nejad et al. (2016) and verified with decay tests. The mean and standard deviation of the fitted error distribution are shown in Fig. 8. The errors in the inertia and stiffness parameter estimation show similar behavior, with maximum errors of σ_SI<17 %. Local maxima in the mean and standard deviation are observed near cut-in (5 m s⁻¹) and near rated wind speeds (11–13 m s⁻¹), while the minimum is located at cut-out wind speed (25 m s⁻¹). It appears that the quasi-stationary conditions in the torque-controlled operational regime above rated wind speeds are conducive to accurate parameter estimation, while the transient dynamics at rated wind speeds due to activation and deactivation of the pitch controller introduce higher estimation errors. Contrary to the inertia and stiffness estimates, the damping parameter estimates show significantly higher errors reaching values of up to σ_SI<105 %. This finding is in agreement with recent studies on drivetrain model validation, where it is reported that the estimation of damping values by OMA techniques is challenging due to the low parameter sensitivity (Vanhollebeke et al., 2015). The damping parameter has a small influence on the dynamic response outside of the resonance area in the considered operational conditions.

Figure 8 Mean and standard deviation of system identification error distributions e_SI.

3.5 Modeling errors

Lastly, the modeling errors e_model in the gear and bearing load calculation due to the complexity reduction of the ROMs are investigated. The discussion is divided into a frequency analysis (Sect. 3.5.1), the analysis of the model bias (Sect. 3.5.2), and the analysis of the dynamic model error (Sect. 3.5.3).

3.5.1 Characterization of drivetrain dynamics

A frequency analysis of the simulated drivetrain loads is conducted to identify which aspects of the drivetrain dynamics the ROMs are able to represent well and which aspects are sources of modeling errors. The drivetrain dynamics can be generally characterized as dynamic responses to a variety of both internal and external excitations. These excitations can be further differentiated into torque and non-torque loads, i.e., lateral forces and bending moments (Table 2). External excitations are mainly the result of aerodynamics and are prevalent at low frequencies. Aerodynamic imbalance is present in healthy conditions due to turbulence, wind shear, the vertical wind profile, and the rotor axis tilt or is caused by faulty yaw and pitch misalignment. This results in periodic load variations in the rotor torque, thrust, and bending moments at the rotor frequency 1P (Mehlan et al., 2023). The tower shadow is also known to induce similar torque and non-torque excitations at the blade passing frequency 3P. The system boundaries of the drivetrain models cut through the rotor hub and the yaw bearing; hence, all structural dynamics of the blades and the tower are considered to be external excitations. These are simulated with the global aeroelastic models, and the resulting main shaft loads and tower motions are applied as boundary conditions in the drivetrain models. The deformation of the blades with edgewise bending modes translates to torque excitations at the main shaft, while flapwise bending modes cause primarily non-torque excitations. Similarly, fore–aft and side–side tower bending introduces excitations in the thrust and bending moments. Internal excitations are caused by periodic changes in component stiffnesses and occur generally at much higher frequencies. Gear mesh excitations are a result of the changing number of tooth contacts during one meshing cycle. Gear meshing primarily results in periodic variation of the transmitted torque but may also have non-torque components in helical gear stages. Bearing excitations are caused by roller elements passing the load zone and result in non-torque excitations at the ball passing frequencies. Further internal excitations are observed at the planet carrier rotational frequencies. Shaft misalignment, mass imbalance, or non-torque loading may result in bending of the flexible planet carrier and in skewing of the load distribution between planets such that each planet bearing experiences periodic load changes during one planet carrier revolution.

Figure 9 Power spectral densities of bearing radial loads simulated with the 5 MW FOM, rigid ROM, and flexible ROM at EC8.

Table 2 Types of excitations and characteristic frequencies in wind turbine drivetrains.

The characteristic excitations are observable in the power spectral densities (PSDs) of the bearing loads (Fig. 9). Shown are the simulated bearing loads at each gear stage for EC8 (17 m s⁻¹) using the FOM and the rigid and flexible ROM. The rigid ROM exhibits good agreement in the lowest frequency range (<1 Hz) governed by wind and wave load excitations but generally underestimates higher-frequency dynamics, as it only considers rigid-body modes. The flexible ROM achieves more accurate load estimates by inclusion of the first torsional drivetrain mode. It is able to match the peaks of external excitations such as the first collective edgewise blade bending mode (f_N1) and higher-order modes. The internal dynamics are captured reasonably well, with good agreement in the second-stage gear meshing frequency (f_gm2). However, some discrepancies remain in the first-stage gear meshing frequency peak (f_gm1) and in the planet carrier excitations (f_plc1, f_plc2), visible at the first- and second-stage planet bearings (PL-A, IMS-PL-A). These suggest the presence of non-torque loads at the planet carriers. The investigated 5 and 10 MW drivetrain models are designed with a four-point main bearing suspension, where it is generally assumed that all non-torque loads of the rotor are fully compensated for by the main bearings, but it appears that this is not the case and that non-torque loads partially propagate further downwind into the drivetrain. The results showcase the limitations of torsional ROMs and suggest that a significant source of modeling errors originates from neglecting planetary carrier bending modes.

Figure 10 Mean value of fitted modeling error distributions e_model indicating the ROM's bias.

3.5.2 Model bias

The focus of the statistical analysis lies first on the model bias quantified by the mean of the fitted error distribution. Shown in Fig. 10 are the model biases of the rigid and flexible ROM in numerical and experimental case studies. The field measurements are only available for the HSS-A bearing. The highest biases are observed near cut-in wind speeds (5 m s⁻¹), which can be associated with start-up and shut-down effects. At higher wind speeds (> 7 m s⁻¹) the environmental conditions have a marginal influence on the model bias. Significant biases of up to 36 % are observed at the high-speed gear stage. The loads at the upwind HSS-A and IMS-A bearings are consistently underestimated, while the loads at the downwind HSS-B and IMS-B bearings are overestimated. One reason for these discrepancies could lie in the physical simplifications of the ROMs, which reduces the gear contact force to a singular vector along the line of action. The load distribution along the gear flank is not considered, and thus the bending moments resulting from inhomogeneous load distributions are neglected. Other authors introduce a “twist stiffness” perpendicular to the circumferential gear meshing stiffness to account for the load distribution (Eritenel and Parker, 2012). However, in this approach the solution requires knowledge of gear and bearing stiffness parameters, which are difficult to determine and validate in practice. Another factor could be the assumption of open-ended shafts that do not allow the transfer of non-torque loads. In the FOMs this is not the case, since the generator coupling at the HSS and the sun–planet gear contact at the IMS allow the transfer of shear forces. These could skew the HSS and IMS bearing loads and further contribute to the model bias. The persistence of model biases in such analytical ROMs is further supported with field measurements of the DOE 1.5 MW turbine. The measured model bias is independent of the EC and amounts to about 12 %, which is of similar magnitude as the values of the numerical case studies.

Figure 11 Standard deviation of fitted modeling error distributions e_model indicating the ROM's dynamic errors.

3.5.3 Dynamic error

The standard deviation of the of the fitted error distributions indicates how well the ROMs capture drivetrain dynamics compared to the FOM. As depicted in Fig. 11, σ_model is positive for all considered cases, which suggests that the ROMs generally underestimate the load dynamics. The error distributions show similar trends across all bearing and gear types. The highest values are observed near cut-in wind speeds (5 m s⁻¹), followed by a steep decline to the global minimum at 9 m s⁻¹ and a gradual progressive trend towards cut-out wind speeds (25 m s⁻¹). Similarly to the high model bias, the high errors at cut-in wind speeds can be attributed to start-up and shut-down effects. The progressive trend can be attributed to aerodynamic non-torque loads transferred from the rotor into the drivetrain. While the torque is controlled to rated conditions above rated wind speed, the non-torque loads, in particular pitch and yaw bending moments, continue to increase with higher wind speeds (Mehlan et al., 2022b). These can excite non-torsional modes of the drivetrain, in particular planet carrier bending modes (see Sect. 3.5.1), which the purely torsional ROMs do not account for. The flexible ROM appears to capture the drivetrain dynamics to a much higher degree than the rigid ROM, resulting in lower modeling errors across all bearing and gear locations. The largest differences are observed above rated wind speed, where the excitation of the first drivetrain torsional mode becomes increasingly more energetic. Below rated wind speed the relative improvement is much lower, since in this operational regime the drivetrain dynamics are governed by rigid-body modes. The modeling errors observed in the field measurements show a similar trend and order of magnitude and support the previous findings of the numerical case studies.

3.6 Long-term fatigue damage error

The use case of long-term fatigue damage monitoring is considered to assess the impact of the modeling and estimation errors in the DT framework. Three scenarios are hereby considered with increasing resolution of SCADA measurements, ranging from 10 min and 1 to 200 Hz. The resolution of 10 min and 1 Hz limits the DT model to the rigid torsional ROM, since the first torsional natural frequency lies above the Nyquist frequency, while the case of 200 Hz measurements allows the application of the flexible ROM. The long-term fatigue damage is calculated by weighting the short-term fatigue damage of each EC with the wind speed distribution. As shown in Fig. 12, the contribution of wind speeds near cut-in (3–7 m s⁻¹) to long-term fatigue does not exceed 2 % due to the low probability of such wind speeds in addition to small aerodynamic loads. The small contribution suggests that the high errors observed at cut-in wind speeds due to start-up and shut-down effects (Sect. 3.5.2) have a negligible impact. The highest contribution has wind speeds of 13 m s⁻¹, where model and measurement errors are fortunately near their minima. The relative error in long-term fatigue damage for each of the scenarios is shown in Fig.13. The long-term fatigue damage is generally underestimated by the DTs due to underestimation of the load amplitudes. It should be noted that the error in the bearing and gear load estimates is amplified by exponentiation with the S–N curve exponent of 10 $/$ 3 and 6.225, respectively. Hence, the gear fatigue damage error tends to be larger due to the larger exponent. The first scenario with 10 min SCADA data results in relative errors of up to −44.4 % in the gear fatigue damage and up to −15.9 % in the bearing fatigue damage due to the high measurement errors e_meas (Sect. 3.2). The resolution is insufficient to capture the low-frequency aerodynamics and the high-frequency internal drivetrain dynamics. The second scenario with 1 Hz data yields significantly smaller relative errors limited to −11.2 % and −6.6 % in the gear and bearing fatigue damage, respectively. In this case, the rigid ROM is able to represent low-frequency load variations due to wind and wave excitations but is limited with respect to higher-frequency internal dynamics. The third scenario with 200 Hz measurements and the two DOF flexible ROM results in only marginally lower fatigue damage errors of −9.7 % and −5.5 %, which showcase the trade-off of increasing the model fidelity. While the addition of a torsional DOF in the flexible ROM significantly reduces the modeling errors e_model (Sect. 3.5.3), it introduces one unknown variable in the rotor torque and four unknown parameters in the rotor inertia, generator inertia, drivetrain stiffness, and damping. The estimation of the rotor torque and the parameters by inverse methods causes additional estimation errors (e_SE, e_SI) (Sect. 3.3 and 3.4), which partially diminish the benefit of the higher model accuracy.

Figure 12 Contribution of each wind speed bin to long-term fatigue damage for the example of the 5 MW HSS-A bearing.

Figure 13 Relative error [%] in long-term bearing and gear fatigue damage.

Table 3 Summary of the error distributions of different DT components for fatigue damage calculation.

4 Conclusions

This paper presents a systematic assessment of the accuracy of DTs for load and fatigue damage monitoring in wind turbine drivetrains. Numerical studies with the NREL 5 MW and DTU 10 MW reference turbines and experimental studies with the DOE 1.5 MW research turbine were conducted to assess modeling and estimation errors of different DT components and their impact on long-term fatigue damage. The information loss in the SCADA data input e_meas, the errors of the state estimation and system identification methods (e_SE, e_SI), and the modeling errors of the drivetrain ROMs e_model were investigated and quantified using normal distributions (Table 3). The investigation of the measurement errors revealed a significant loss of information by using 10 min averaged SCADA data. The measurement resolution is insufficient to observe the low-frequency drivetrain load dynamics due to the wind speed and rotor torque volatility, which resulted in maximum errors of σ=74 % and long-term fatigue damage errors of up to −44.4 % in the gears and −15.9 % in the bearings. The results strongly advocate for the use of high-frequency SCADA data with a resolution of at least 1 Hz for fatigue monitoring purposes. The second source of errors is identified in the state estimation method, the augmented Kalman filter, that is applied to match the dynamic state of the DT model with the physical wind turbine based on real-time data streams. The challenge lies in estimating the rotor torque, which is not measured directly and must be estimated by the Kalman filter. The Kalman filter tends to underestimate the rotor torque at the first torsional natural frequency, which results in errors ranging σ=1 %…7 % at normal operational conditions (> 5 m s⁻¹). The third source of errors originates from the aleatory uncertainty of the system properties. Inertia, stiffness, and damping values may vary over the turbine's life cycle as a result of faults, material degradation, or part replacement. System identification methods are applied to detect these changes and update the model parameter accordingly. The errors in the parameter estimates are particularly high at cut-in and near rated wind speeds (σ<17 %) due to transient dynamics and the high variance in the drivetrain torque, while the lowest estimation errors are observed in the torque-controlled regime above rated wind speed (σ>0.5 %). Furthermore, it is observed that the estimation of the drivetrain torsional damping is significantly more inaccurate than inertia and stiffness parameters (σ<105 %). This is likely due to the low sensitivity of the damping parameter with respect to the drivetrain torsional dynamics at normal power production. Lastly, the modeling errors due to the ROMs' limitations is investigated. ROMs with one or two torsional DOFs are used as DT models due to their lower computational loads in real-time monitoring applications, their lower validation costs, and the limited observability of non-torsional dynamic states with the available SCADA measurements. One-DOF rigid ROMs are only able to match the dynamics in the lowest frequency range (<1 Hz) governed by wind and wave load excitations, while two-DOF flexible ROMs better capture the dynamic drivetrain response to higher-frequency internal excitations such as gear meshing. Remaining limitations are observed in capturing non-torsional dynamics, in particular the bending dynamics of the first- and second-stage planet carriers. The load estimation errors of the flexible ROM are noticeably smaller; however, only a small improvement with respect to the fatigue damage estimates is observed (−6.6 % to −5.5 %). While the addition of a torsional DOF in the flexible ROM significantly reduces the modeling errors, it introduces additional unknown variables and parameters with associated estimation errors that partially diminish the benefit of the higher model accuracy. The presented study contributes to a deeper understanding of the modeling and estimation errors in DTs for load and fatigue monitoring. The error distributions may be used in reliability studies, in risk assessment and the derivation of safety factors, or in decision processes on model fidelity and sensor measurement resolution.