Identification of wind turbine main-shaft torsional loads from high-frequency SCADA (supervisory control and data acquisition) measurements using an inverse-problem approach

To assess the structural health and remaining useful life of wind turbines within wind farms, the site-specific structural response and modal parameters of the primary structures are required. In this regard, a novel inverse-problem-based methodology is proposed here to identify the dynamic quantities of the drivetrain main shaft, i.e. torsional displacement and coupled stiffness. As a model-based approach, an inverse problem of a mathematical model concerning the coupled-shaft torsional dynamics with high-frequency SCADA (supervisory control and data acquisition) measurements as input is solved. It involves Tikhonov regularisation to minimise the measurement noise and irregularities on the shaft torsional displacement obtained from measured rotor and generator speed. Subsequently, the regularised torsional displacement along with necessary SCADA measurements is used as an input to the mathematical model, and a model-based system identification method called the collage method is employed to estimate the coupled torsional stiffness. It is also demonstrated that the estimated shaft torsional displacement and coupled stiffness can be used to identify the site-specific main-shaft torsional loads. It is shown that the torsional loads estimated by the proposed methodology is in good agreement with the aeroelastic simulations of the Vestas V52 wind turbine. Upon successful verification, the proposed methodology is applied to the V52 turbine to identify the site-specific main-shaft torsional loads and damage-equivalent load. Since the proposed methodology does not require a design basis or additional measurement sensors, it can be directly applied to wind turbines within a wind farm that possess high-frequency SCADA measurements.

1 Introduction

Monitoring of wind turbines within wind farms is increasingly becoming very important due to the need to detect anomalous behaviour, plan inspections or preventive maintenance, and compute the remaining useful life of specific structures. The site-specific structural dynamic quantities such as structural response and modal parameters could assist in the condition monitoring of wind turbines.

The structural response of the wind turbines is measured using load instrumentation such as accelerometers (Koukoura et al., 2015; Pahn et al., 2017; Norén-Cosgriff and Kaynia, 2021), and an output-based operational modal analysis (OMA) (Wang et al., 2016) technique is employed for the identification of the modal parameters. OMA methods can be broadly classified into two categories: (i) time-domain-based methods and (ii) frequency-domain-based methods (Zhang et al., 2010). A recent comprehensive review on various time-domain- and frequency-domain-based OMA methods can be found in Zahid et al. (2020). Upon identification of the modal parameters together with the measured structural response, an inverse-problem-based technique is employed for the estimation of the site-specific loads (Pahn et al., 2017). Alternatively, one can use strain-gauge-based load sensors to measure the site-specific loads directly (IEC 61400-13, 2016).

The addition of new instrumentation to existing turbines, such as the installation of strain gauges and accelerometers, can be costly and also require repetitive calibration and synchronisation of their measurement signals with the turbine computer. Most wind farm operators also do not possess the aeroelastic design parameters of their wind turbines and hence cannot simulate the mechanical loads acting on their wind turbines. Monitoring turbine primary structures through existing SCADA-based (supervisory control and data acquisition) measurements can yield a cost-effective solution and provide valuable information to the wind farm operator. Usually, such monitoring through SCADA only provides information on power performance without regarding turbine structural integrity. Since there are two SCADA signals (i.e. rotor speed and generator speed) related to torsional oscillations of the main shaft of wind turbine drivetrains, the same can be used to quantify the torsional dynamics of the main shaft.

For this purpose, an inverse-problem-based approach is developed here to determine the torsional stiffness and response of the main shaft of a wind turbine, using existing high-frequency SCADA measurements such as the rotor speed and generator speed. This is a cost-effective alternative approach that is being proposed for the main shaft without using any additional measurement sensors or an aeroelastic design basis of the wind turbine. The proposed inverse-problem approach is a model-based approach whereby a mathematical model concerning the shaft torsional dynamics will be utilised to obtain both the torsional displacement and the coupled torsional stiffness in a continuous time domain. It involves Tikhonov regularisation (Tikhonov, 1963) for regularising the measurement data and the collage method (Kunze and Vrscay, 1999) for estimating the torsional stiffness.

The concerned mathematical model comprised of differential equations will be solved for the shaft torsional velocity with high-frequency rotor and generator speed measurements as inputs. Subsequently, the main-shaft torsional displacement is obtained by numerically integrating the shaft torsional velocity. However, numerical integration is based on time-marching algorithms, and the lack of initial conditions makes displacement reconstruction an ill-posed problem (Hong et al., 2008). Since it is an ill-posed problem, influence of measurement noise will be amplified during the time-marching procedure, resulting in an erroneous displacement. This inaccurate displacement also leads to drastic errors in the system modal-parameter estimation. Hence, one needs to go for regularisation techniques to smoothen the reconstructed displacement. Hansen (2005) discussed the nature of various ill-posed problems and presented several solution methodologies. Though there are many regularisation techniques available, Tikhonov, truncated singular value decomposition (SVD), and nuclear norm are a few of the popular techniques (Aarden, 2017). Among all these regularisation techniques, Tikhonov regularisation has been widely used in many engineering applications (Ronasi et al., 2011; Hào and Quyen, 2012; Bangji et al., 2017; Nieminen et al., 2011), and it has been studied extensively in the field of inverse problems (Hansen, 2005) as well. Further, digital filters and the frequency domain integration approach (FDIA) are also widely used techniques in the literature to reconstruct displacements from measured accelerations (Hong et al., 2008; Brandt and Brincker, 2014; Qihe, 2019; Lee et al., 2010). However, digital filters such as infinite-impulse-response (IIR) and finite-impulse-response (FIR) filters have several drawbacks when reconstructing the low-frequency dominant displacements as is the case here (Lee et al., 2010). On the other hand, the FDIA methods are sensitive to the time interval of the measurements (Lee et al., 2010). It is shown by Lee et al. (2010) that Tikhonov regularisation is better suited for low-frequency dominant structures. Hence, the same has been employed in the present work. Tikhonov regularisation minimises the error using the least-square criterion and by means of numerical damping; it also minimises the effect of measurement noise.

Upon obtaining the regularised shaft torsional displacement, the same mathematical model is utilised to obtain the shaft stiffness. For this purpose, the collage method – a model-based system identification technique – is used (Kunze and Vrscay, 1999; Groetsch, 1993). This method has been successfully applied for the system identification in various differential-equation-based problems such as boundary value problems (Kunze et al., 2009), reaction–diffusion problems (Deng et al., 2008), and elliptic problems (Kunze and La Torre, 2016). The collage method transforms the system identification problem into a minimisation problem of a function of several variables (for example unknown system parameters), and then the minimisation problem is solved using a minimisation algorithm called collage coding (Kunze et al., 2009). Collage coding is a greedy algorithm that attempts to find the approximate solution in a single step without any need for iteration, as is the case for other inverse-problem-based methods (Deng and Liao, 2009). Further, the model-based collage method is simple, easy to implement, and computationally inexpensive as compared to the output-based OMA methods.

The estimated shaft torsional stiffness and displacement are further used to identify the site-specific shaft torsional load. This novel methodology can potentially benefit wind farm owners, since both the property of the structure in terms of its stiffness and the structural response and the site-specific load can be determined without requiring additional sensors or information from the wind turbine manufacturer. The main-shaft torsional load affects the fatigue performance of other drivetrain components such as the gearbox and planetary bearings (Dong et al., 2012; Gallego-Calderon and Natarajan, 2015; Ding et al., 2018). Hence, the same site-specific torsional load may be used as an input for quantifying the remaining useful life (RUL) (Ziegler et al., 2018) of the main shaft, the gearbox, and other drivetrain components as well. The estimation of RUL and yearly damage does not require additional historical weather data and condition monitoring data, as the wind speed and wind direction measurements are available in the SCADA measurements. However, this is beyond the scope of present work. Further, the proposed approach requires that the sampling frequency of the SCADA measurement be significantly higher than the dominant frequencies of the drivetrain torsional oscillations (i.e. 1p and 3p rotor excitation frequencies and torsional natural frequencies, where p is the mean rotor speed). As a result, the proposed method cannot be used for the turbines that have measurements in terms of 10 min SCADA statistics.

The rest of the paper is organised as follows: the problem formulation consisting of Tikhonov regularisation and the collage method is given in Sect. 2; Sect. 3 presents the verification of the proposed formulation; and application of the proposed formulation on measurements is presented in Sect. 4.

Figure 1 A two mass model of the wind turbine drivetrain.

2 Problem formulation

As mentioned in the previous section, the main objective is to identify the shaft torsional displacement and coupled stiffness from SCADA measurements. This is achieved by solving the shaft torsional dynamical equations using a suitable inverse-problem algorithm and the estimated shaft torsional displacement θ, and torsional stiffness K will be utilised for the shaft torsional-load estimation. For this purpose, a two-mass model (refer to Fig. 1) (Boukhezzar et al., 2007; Girsang et al., 2013; Berglind et al., 2015) that governs the main-shaft torsional dynamics subjected to the rotor and generator torques T_r and T_g, respectively, is considered, and the mathematical model is given by Eqs. (1)–(3). The governing equations are obtained in an inertial frame of reference and converted to the low-speed side of the drivetrain by means of the gear ratio (Boukhezzar et al., 2007). By assuming that the drivetrain components are in a series representation in terms of its modal quantities, effective values for the modal quantities are used in the two-mass model (Girsang et al., 2013).

$$\begin{array}{}\text{(1)}& {\displaystyle }{J}_{\mathrm{r}}{\dot{\mathit{\omega }}}_{\mathrm{r}}={\mathit{T}}_{\mathrm{r}}-K\mathit{\theta }-C\dot{\mathit{\theta }},\text{(2)}& {\displaystyle }{J}_{\mathrm{g}}{\dot{\mathit{\omega }}}_{\mathrm{g}}=-{\mathit{T}}_{\mathrm{g}}+K\mathit{\theta }+C\dot{\mathit{\theta }},\text{(3)}& {\displaystyle }\dot{\mathit{\theta }}={\mathit{\omega }}_{\mathrm{r}}-{\mathit{\omega }}_{\mathrm{g}}.\end{array}$$

Here, J_r represents the inertia of the rotor; J_g represents the collective inertias of the high-speed shaft (HSS), the gearbox, and the generator; ω_r and ω_g are the rotor and generator speeds, respectively; K and C are the effective stiffness and damping of the drivetrain including the main shaft, HSS, and gearbox; and θ is the torsional displacement of the drivetrain. Throughout the article, all vector quantities have been marked in bold.

Given the modal parameters ($J_{\mathrm{r}},J_{\mathrm{g}},C$, and K) and external torques (T_r and T_g), Eqs. (1)–(3) are solved for ω_r,ω_g, and θ, which is known as a forward-problem approach (Pahn, 2013). But given only SCADA measurements, Eqs. (1)–(3) are solved inversely for θ and modal parameters. The available SCADA measurements are ${\mathit{\omega }}_{\mathrm{r}},{\mathit{\omega }}_{\mathrm{g}},P,\mathit{\beta }$, and U. Here, P,β, and U are, respectively, the generator power, the blade pitch angle, and wind velocity. The proposed inverse-problem approach consists of Tikhonov regularisation for regularising the measurement data and the collage method for estimating the torsional stiffness, and the entire methodology is shown as a flowchart in Fig. 2.

Figure 2 Flowchart depicting the inverse-problem algorithm.

In the following, implementation of Tikhonov regularisation and the collage method on the drivetrain torsional dynamics will be discussed in detail.

2.1 Tikhonov regularisation

Given ω_r and ω_g, $\dot{\mathit{\theta }}$ is obtained by using Eq. (3), and then by numerically integrating $\dot{\mathit{\theta }}$, the shaft torsional displacement (θ) is obtained. The numerical-integration schemes require initial conditions as they march on time. However, the initial conditions are unavailable or inaccurate in practice. The lack of initial conditions (usually assumed to be 0) on the displacement are inconsistent with the real values that result in the phenomenon of baseline shift or drift which causes the position error to grow with time during the integration (Pahn et al., 2017). The effect of the lack of initial conditions on the reconstructed displacement obtained using the time integration scheme is shown in Fig. 3. As seen in the figure, the numerical error is multiplicatively increased with time, which results in a drift in the reconstructed displacement. As mentioned in the introduction, to minimise the numerical error due to the lack of initial conditions and to minimise the effect of measurement noise, a widely used regularisation technique called Tikhonov regularisation (Tikhonov, 1963) is employed here.

Figure 3 Comparison of time integration displacement with actual displacement. The torsional displacement is in radians.

Implementation of Tikhonov regularisation on the velocity to obtain the displacement is not readily available in the literature, and hence the same is presented in Appendix A1 for the sake of completeness. By following the procedure outlined in Appendix A1, the regularised torsional displacement (θ) is obtained as

$$\begin{array}{}\text{(4)}& \mathit{\theta }=({\displaystyle \frac{{\mathbf{L}}^{\mathrm{2}}}{\mathrm{4}}}+{\mathit{\lambda }}^{\mathrm{2}}\mathbf{I}{)}^{-\mathrm{1}}{\displaystyle \frac{{\mathbf{LL}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t}{\mathrm{2}}}.\end{array}$$

Here, the matrices L,L_a, and I are defined in the Appendix A1; λ is the regularisation parameter; Δt is the time interval; and $\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}$ is the torsional velocity obtained from Eq. (3). The obtained regularised torsional displacement is compared with the actual displacement in Fig. 4 along with the numerical-integration result. As seen in Fig. 4, there is a close match between the result by Tikhonov regularisation and the actual displacement as compared to the numerical-integration result.

Figure 4 Comparison of the Tikhonov and time integration displacements with actual displacement.

2.2 Collage method

Upon estimating θ using Tikhonov regularisation, the next step is to estimate the modal parameters using the collage method (Kunze and Vrscay, 1999; Groetsch, 1993) – a model-based approach. The mathematical formulation of the collage method and its implementation are given in Appendix A2. The rotor equation (Eq. 1) cannot be used for the parameter estimation, as there is no information about the rotor torque in the SCADA measurement. Instead, the collage method is employed on the generator equation (Eq. 2) for estimating the modal parameters, as the generator torque (T_g) can be readily obtained from the SCADA data as ${\mathit{T}}_{\mathrm{g}}=P/{\mathit{\omega }}_{\mathrm{g}}$. Accordingly, for the target function ${\mathit{\omega }}_{\mathrm{g}}\left(t\right),t\in [t_{\mathrm{0}},t_n]$, the squared ℒ² collage distance (refer to Appendix A2) for Eq. (2) becomes

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{\Delta }}^{\mathrm{2}}& =\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\left[J_{\mathrm{g}}\right[{\mathit{\omega }}_g-{\mathit{\omega }}_{g_{\mathrm{0}}}]+\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}{\mathit{T}}_{\mathrm{g}}\mathrm{d}t\\ & -K\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\mathit{\theta }\mathrm{d}t-C\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\dot{\mathit{\theta }}\mathrm{d}t{]}^{\mathrm{2}}\mathrm{d}t,\end{array}\end{array}$$

where ${\mathit{\omega }}_{{\mathrm{g}}_{\mathrm{0}}}$ is the generator speed at time, t=t₀. Modal parameters will be obtained by minimising Eq. (5) with respect to J_g,C, and K. Upon obtaining θ using Tikhonov regularisation and K using the collage method, the torsional load is obtained as M_z=Kθ.

2.2.1 Application of the collage method

To test the applicability and efficiency of the collage method for the wind turbine drivetrain system, a verification study is undertaken by comparing the main-shaft torsional stiffness obtained using the collage method with its design value for two different wind turbines. This is done for the following two wind turbines, (i) DTU (Technical University of Denmark) 10 MW (Bak et al., 2013) and (ii) Vestas-V52 (Vestas, 2020). For facilitating a comparison of the main-shaft torsional response alone, the rigid variant of the turbine is chosen, and this implies that the rotor and tower are rigid and that the main shaft alone is considered to be flexible. By performing aeroelastic simulations on these turbines, the shaft torsional displacement θ is obtained. Throughout the article, the aeroelastic simulation is performed using the DTU in-house tool called HAWC2 (Larsen and Hansen, 2007). HAWC2 (Horizontal Axis Wind turbine simulation Code 2nd generation) is used for calculating the wind turbine aeroelastic loads and responses in the time domain. It uses multibody formulation to model the structure, models based on blade element momentum (BEM) theory for modelling the wind effects, and hydrodynamic models for modelling the hydro-effects (in the case of offshore turbines) on the structure. Control of the turbine is performed through the dynamic link libraries (DLLs). Using the main-shaft torsional-load time series obtained from HAWC2 and by minimising Eq. (5) concerning the modal parameters, K,C, and J_g are obtained. Since for the estimation of shaft torsional load, only K is needed among all the modal parameters, the same is compared with the design values. The estimated shaft stiffness and the stiffness from the aeroelastic model of the DTU 10 MW turbine along with percentage error are tabulated in Table 1. Only the percentage error is given for the Vestas V52 turbine in Table 1. As seen in the table, the estimated torsional-stiffness values match well with the design values. If the torsional displacement (θ) is known, then the determination of torsional stiffness (K) from Eq. (5) is readily feasible as explained. However in practice, the shaft torsional displacement is unknown, and therefore the collage equations may not be directly used to determine the shaft stiffness. In the following, the entire proposed methodology will be verified with the aeroelastic simulation results of the Vestas V52 turbine.

Table 1 Comparison of estimated K with design K for a turbine variant where only the main shaft is flexible.

3 Verification of the method on V52 turbine simulations

Aeroelastic simulations are performed for the Vestas V52 turbine corresponding to the design load case (DLC 1.2) (IEC 61400-1, 2019) in HAWC2. The DLC 1.2 was run with 18 10 min load simulations (three yaw directions: 0^∘ ± 10^∘, and six turbulent wind seeds) for each mean wind speed ranging from 4 to 26 m/s in the interval of 2 m/s, which results in a total of 216 simulations. Each simulation uses 10 min normal wind turbulence inflow with a sampling frequency of 50 Hz. From these 216 simulations, the inputs (ω_r,ω_g, and T_g) for the proposed methods are obtained, and by following the procedure depicted in Fig. 2, the torsional loads are obtained, and the same is compared with the simulation results. From the simulated ω_r and ω_g, the torsional velocity $\dot{\mathit{\theta }}$ is obtained using Eq. (3), and then the corresponding θ for each simulation is obtained by using Tikhonov regularisation (Eq. 4). By assuming that the same level of numerical noise (noise due to a lack of initial conditions) presents in all the simulations of a particular mean wind speed, the optimal λ parameter for regularisation (refer to Appendix A1) is estimated only once per mean wind speed. The estimated regularisation parameter (λ) for each mean wind speed is presented in Fig. 5. As seen in the figure, a higher value of λ at higher wind speeds indicates that more numerical damping is needed to suppress the numerical noise.

Figure 5 Estimated λ for each mean wind speed.

Figure 6 Comparison of regularised θ_dyn with the results from aeroelastic simulation. WSP: mean wind speed.

Since the displacement is reconstructed from the velocity, a dynamic quantity, the reconstructed displacement will have a mean of zero, and this displacement component is referred to as a dynamic component of the displacement (θ_dyn). However, the torsional displacement will have a contribution from the static load, and the displacement due to the static load is referred to as a static component of the displacement (θ_stat) (i.e. the mean value of the displacement). The dynamic component of the displacement oscillates about the static component of the displacement. The regularised θ_dyn for two representative mean wind speeds is compared with the dynamic component of simulated torsional displacement in Fig. 6. By removing the static displacement from the simulated displacement, the resulting dynamic component can be compared with the regularised θ_dyn as shown in Fig. 6. As seen in the figure, the regularised θ_dyn matches well with the simulated dynamic displacement for most of the time except for the peak amplitudes.

Upon ensuring the correctness of θ_dyn, the collage method can be employed for K estimation using the results of all 216 simulations. Since only θ_dyn is available, Eq. (5) is modified to have the dynamic components alone for the estimation of K as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{\Delta }}^{\mathrm{2}}& =\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\left[J_{\mathrm{g}}\right[{\mathit{\omega }}_{\mathrm{g}}-{\mathit{\omega }}_{{\mathrm{g}}_{\mathrm{0}}}]+\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}({\mathit{T}}_{\mathrm{g}}-{\mathit{T}}_{{\mathrm{g}}_{\mathrm{mean}}})\mathrm{d}t\\ & -K\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}{\mathit{\theta }}_{\mathrm{dyn}}\mathrm{d}t-C\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\dot{\mathit{\theta }}\mathrm{d}t{]}^{\mathrm{2}}\mathrm{d}t.\end{array}\end{array}$$

By minimising Eq. (6), K will be obtained for all 216 simulations. Owing to uncertainty associated with θ_dyn estimation along with the collage method, three skewed distributions for K can be determined at different mean wind speeds. For these distributions, the mode is a more stable estimate of the central tendency, as it is least biased by outliers (Hedges and Shah, 2003). The mean of modes of the resulting probability density functions (PDFs) provides the resultant estimate of K. The relative error between the estimated K value obtained by taking the mean of modes and the design value of the turbine is 12.06 %. At this point, it is important to remember that since the inputs are from the HAWC2 aeroelastic simulation which does not account for gearbox dynamics, the estimated stiffness value only has the contributions from the rotor and the main shaft. As a result, the estimated K is the resultant torsional stiffness of the main shaft including the blade edgewise stiffness contribution.

Even though the dynamic component of the displacement is sufficient enough for K estimation as explained above, the mean value (i.e. the static component) of the displacement has a significant effect on the fatigue damage (Veldkamp, 2006). Hence, it is important to estimate the same, and this is obtained by solving the static problem of the drivetrain. Considering the static equilibrium of Eq. (1) with all parameters expressed on the low-speed side is as follows:

$$\begin{array}{}\text{(7)}& K{\mathit{\theta }}_{\mathrm{stat}}={\mathit{T}}_{{\mathrm{r}}_{\mathrm{mean}}}={\displaystyle \frac{{\mathit{T}}_{{\mathrm{g}}_{\mathrm{mean}}}}{{\mathit{\eta }}_{\mathrm{gen}}}}.\end{array}$$

Here, η_gen is the generator efficiency, which is 94.4 % for the V52 turbine. As an alternative, one can use the overall efficiency of the wind turbine as well for the sake of completeness. However, the use of different efficiencies will not significantly affect the outcome. From Eq. (7), θ_stat is obtained as

$$\begin{array}{}\text{(8)}& {\mathit{\theta }}_{\mathrm{stat}}={\displaystyle \frac{{\mathit{T}}_{{\mathrm{g}}_{\mathrm{mean}}}}{K{\mathit{\eta }}_{\mathrm{gen}}}}.\end{array}$$

Finally, the static and dynamic components are superimposed to get the actual displacement (θ). This actual displacement will then be used to estimate the site-specific main-shaft torsional load as shown by Eq. (9).

$$\begin{array}{}\text{(9)}& {\mathit{M}}_z=K\mathit{\theta }=K({\mathit{\theta }}_{\mathrm{stat}}+{\mathit{\theta }}_{\mathrm{dyn}})\end{array}$$

At this stage, the estimated torsional loads can be compared with the simulated torsional loads. The comparison of the torsional loads for two representative mean wind speeds is shown in Fig. (7). As seen in Fig. (7), all the important aspects of the time-series variation and the dominant frequency dynamics (low-frequency components – up to first three peaks) are captured quite well in the estimated torsional load. The computational time for identifying the regularisation parameter for each mean wind speed is about 40 min in real time, and the stiffness estimation for the 10 min simulation takes 14 s in real time. The above computations are performed on a single node of the high-performance computing cluster of DTU. The cluster has 320 nodes in total, with each node consisting of two Intel Xeon E5-2680 v2 processors, and each processor consists of 10 cores running at 2.8 GHz. If the wind farm is of same type of turbine, then it is sufficient to estimate the stiffness for one of the turbines. However, the shaft torsional displacement needs to be estimated for all the turbines.

Figure 7 Comparison of reconstructed time series and power spectral density (PSD) of the torsional load for two different mean wind speeds.

Upon estimating the torsional load, the torsional damage-equivalent load (DEL) at each mean wind speed is calculated using the following equation:

$$\begin{array}{}\text{(10)}& \text{DEL}=\left({\displaystyle \frac{\mathrm{1}}{N_{\mathrm{ref}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{sim}}}\right({\displaystyle \frac{\mathrm{1}}{N_{\mathrm{sim}}}})\sum _{k=\mathrm{1}}^{k_n}{\displaystyle \frac{N_{i,k}{S}_{i,k}^m\left(\mathrm{0}\right)}{T_{\mathrm{sim}}}}{)}^{\frac{\mathrm{1}}{m}},\end{array}$$

where T_sim is the duration of the load case, N_sim is the number of simulations at each mean wind speed, k_n is the total number of load cycles in a given time series, N_i,k are the number of cycles at load amplitude S_i,k(σ) determined with rain flow counting, and m is the Wöhler exponent. The zero-mean load amplitude is obtained as $S_{i,k}\left(\mathrm{0}\right)=S_{i,k}\left(\mathit{\sigma }\right)+M{S}_{\mathrm{m}}$ (Veldkamp, 2006), where S_m is the mean load and M is the mean stress sensitivity. The turbine shaft is made up of cast iron, and Veldkamp (2006) has reported that for such material, the mean stress correction factor is M=0.19. The Wöhler exponent for the cast iron of m=6 (Veldkamp, 2006) is used in the present study. When N_ref in Eq. (10) equals the component's lifetime in seconds, the DEL has the frequency of 1 Hz. The computed torsional DEL using Eq. (10) for DLC 1.2 is compared with the simulated DEL and is shown in Fig. (8). For most of the mean wind speeds, the estimated DEL is in good agreement with the simulated DEL, and the absolute error between these two at each mean wind speed ranges from 4 % to 12 %. Higher error at higher mean wind speeds is due to the fact that the peak amplitudes in θ_dyn are not captured well using Tikhonov regularisation. The fact that the peak amplitudes are not captured in the reconstructed torsional displacements is typical for Tikhonov regularisation, as there is a slight mismatch in the frequency spectra between the actual and reconstructed torsional oscillations (Lee et al., 2010). Another reason could be due to the fact that the pitch angle influences the main-shaft torsional oscillation through the rotor torque and the rotor speed. However, the controller is maintaining a constant rotational speed beyond the rated wind speed; hence the instantaneous changes in the pitch angle are not propagated through the rotor speed to the main-shaft oscillations. All these effects, in addition to the uncertainty in K estimation, may have resulted in a maximum error of 12 % in the estimated DEL.

Figure 8 Comparison of the predicted DEL with the DEL computed from aeroelastic simulations over all mean wind speeds.

Figure 9 Regularised θ_dyn for two mean wind speeds.

Figure 10 Identified torsional loads (M_z) for two mean wind speeds.

Figure 11 Estimated torsional DEL from SCADA measurements of the Vestas V52 turbine.

4 Application on V52 turbine measurements

Upon verifying the proposed method, the drivetrain main-shaft torsional loads are estimated from the SCADA measurements without any need for the aeroelastic model. SCADA measurements taken during January 2019 for the Vestas V52 850 kW research turbine installed at the DTU Risø site is utilised for this purpose. The measurement data consist of 4459 10 min recorded cases with 50 Hz sampling frequency. Generator torque is used as one of the SCADA signals instead of the generator speed for this part of the study, and the generator speed is obtained from the generator power and generator torque (on the low-speed side) SCADA signals as ${\mathit{\omega }}_{\mathrm{g}}=P/{\mathit{T}}_{\mathrm{g}}$. Further, the measurement data correspond to normal operation and are filtered based on the conditions given in Table 2 that result in 627 cases. It is important to note that based on the site sector assessment, the V52 turbine is in wake-free condition between 280 and 320^∘ wind directions.

Table 2 Normal operation filter conditions.

Using the rotor speed and generator speed, θ_dyn is calculated using the Tikhonov regularisation method, and the obtained θ_dyn for two representative mean wind speeds are shown in Fig. 9.

The regularisation parameter is obtained for each mean wind speed measurement using the L-curve criterion and the obtained values not presented here for the sake of brevity. Subsequently, by applying the collage method on Eq. (6) for all 627 cases, the K values are estimated for each load case, and then the resultant K value is obtained by taking the mean of modes of the resulting PDF as described in the previous section. At this point, it is important to remember that the inputs are from measurement; hence the estimated K is a collective stiffness that has a contribution from all drivetrain components including the gearbox. With the estimated K value, θ_stat is calculated for each load case using Eq. (8), and then the torsional load is obtained from Eq. (9). The calculated torsional loads for two representative wind speeds are shown in Fig. 10a and b.

After estimating the torsional loads for all 627 cases, the identified loads are grouped according to the mean wind speeds ranging from 6 to 22 m/s which are subdivided into nine wind speed bins of 2 m/s width each. Subsequently, the torsional DEL is calculated for each mean wind speed using Eq. (10), and the same is shown in Fig. 11. It is important to note that the DEL given in Fig. 8 (simulation) is the design load and that the DEL given in Fig. 11 is the site-specific load. The difference between these two can give an estimate about the remaining torsional fatigue life of the main shaft under normal operating conditions (Ziegler et al., 2018). Further, the estimated main-shaft torsional DEL may be used to quantify the RUL of the drivetrain (Pagitsch et al., 2020). However, this is beyond the scope of the current work. It is also feasible to reconstruct the torsional-load time series, which may be used an input to gearbox design tools to predict the loading within the gearbox.

5 Summary

A novel inverse-problem-based approach is developed for estimating the main-shaft torsional displacement and stiffness by using high-frequency SCADA measurements. A mathematical model describing the coupled-shaft torsional dynamics is used for this purpose. Numerical errors and the effect of measurement noise on the torsional displacement reconstruction are minimised through the Tikhonov regularisation technique. Subsequently, the collage method is used to estimate the main-shaft coupled torsional stiffness. The estimated main-shaft quantities are then used to identify the main-shaft site-specific torsional load. The proposed formulation is successfully verified for the main-shaft torsional loads with the aeroelastic simulation of the Vestas V52 turbine. Upon verification, the methodology is extended to identify the site-specific main-shaft torsional loads of the same turbine by using SCADA measurements. For this purpose, the measurement data from the DTU Risø site are utilised, and the measurement data are filtered and calibrated for the turbine normal operation. Using the identified torsional loads, the torsional DEL is obtained. Depending on the prior information about the stiffness value, one can either use the entire proposed methodology or follow the torsional displacement estimation part of the proposed methodology for the torsional-load identification. Since the site-specific SCADA measurements are used in the analysis, the obtained loads can give a better estimate of the remaining fatigue life of drivetrain components. Monitoring the estimated loads can help in inspection planning and scheduling maintenance activities. As the proposed methodology does not require any design basis or an aeroelastic design basis, it can be used for wind turbines that possess high-frequency SCADA measurements for the estimation of the main-shaft torsional load and DEL.

Appendix A

A1 Displacement reconstruction using Tikhonov regularisation

By definition, the velocity $\dot{\mathit{\theta }}$ is expressed as

$$\begin{array}{}\text{(A1)}& \dot{\mathit{\theta }}\left(t\right)={\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}\approx \dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\left(t\right),\end{array}$$

where $\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\left(t\right)$ is the velocity obtained from Eq. (3) which can be considered a measured velocity. As explained in Sect. 2.1, the lack of initial conditions in addition to the measurement noise leads to erroneous displacement. In order to minimise the error between the actual and measured velocities in the least-square sense, the following minimisation problem has to be solved:

$$\begin{array}{}\text{(A2)}& \text{Min}{\mathrm{\Pi }}_E\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\left(\dot{\mathit{\theta }}\right(t)-\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}(t){)}^{\mathrm{2}}\mathrm{d}t.\end{array}$$

Here, $\dot{\mathit{\theta }}$ is the calculated velocity. In other words, Eq. (A2) gives a measure of how well the actual velocity approximates the measured velocity. By means of the trapezoidal rule, Eq. (A2) is discretised as follows (Hong et al., 2008):

$$\begin{array}{}\text{(A3)}& {\mathrm{\Pi }}_E\left(\mathit{\theta }\right)\approx \Vert {\mathbf{L}}_{\mathrm{a}}(\dot{\mathit{\theta }}-\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}){\Vert }_{\mathrm{2}}^{\mathrm{2}}\mathrm{\Delta }t,\end{array}$$

where Δt is the time interval of the discretisation and L_a is the diagonal weighing matrix of order (N+1) as

$$\begin{array}{}\text{(A4)}& {\mathbf{L}}_{\mathrm{a}}=\left[\begin{array}{ccccc}\mathrm{1}/\sqrt{\mathrm{2}}& & & & \\ & \mathrm{1}& & & \\ & & \mathrm{\ddots }& & \\ & & & \mathrm{1}& \\ & & & & \mathrm{1}/\sqrt{\mathrm{2}}\end{array}\right].\end{array}$$

Here, N is the number of data points in θ, and for a 10 min simulation with 50 Hz sampling frequency N becomes 30 000. Further, the calculated velocity $\dot{\mathit{\theta }}$ is discretised by the central-difference rule and written in matrix form as

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\mathbf{L}}_{\mathrm{c}}\mathit{\theta }=\dot{\mathit{\theta }},\end{array}$$

where the central-difference matrix L_c of size $(N+\mathrm{1})\times (N+\mathrm{3})$ and the displacement vector θ of size (N+3) are given as

$$\begin{array}{}\text{(A6)}& {\mathbf{L}}_{\mathrm{c}}=\left[\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{\dots }& \\ & -\mathrm{1}& \mathrm{0}& \mathrm{1}& \\ & & \mathrm{\ddots }& & \\ & -\mathrm{1}& \mathrm{0}& \mathrm{1}& \\ & \mathrm{\dots }& -\mathrm{1}& \mathrm{0}& \mathrm{1}\end{array}\right],\mathit{\theta }=\left(\begin{array}{c}{\mathit{\theta }}_{-\mathrm{1}}\\ {\mathit{\theta }}_{\mathrm{0}}\\ \mathrm{⋮}\\ {\mathit{\theta }}_n\\ {\mathit{\theta }}_{n+\mathrm{1}}\end{array}\right).\end{array}$$

In the finite-difference discretisation of Eq. (A6), some defined nodes are located outside of the domain considered (i.e. domain here is time, $t\in [t_{\mathrm{0}},t_n]$ satisfying $t_{\mathrm{0}}\le t\le t_n$). These nodes defined by time steps $i=-\mathrm{1}$ and $i=n+\mathrm{1}$ are fictitious. These nodes are dummy nodes that are used in solving the differential equations by the finite-difference method (Lapidus and Pinder, 2011). Substitution of Eq. (A5) into Eq. (A3) leads to

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}\text{Min}{\mathrm{\Pi }}_E\left(\mathit{\theta }\right)& \approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }t}}{\mathbf{L}}_{\mathrm{a}}{\mathbf{L}}_{\mathrm{c}}\mathit{\theta }-{\mathbf{L}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}{\Vert }_{\mathrm{2}}^{\mathrm{2}}\mathrm{\Delta }t\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{L}\mathit{\theta }-{\mathbf{L}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t{\Vert }_{\mathrm{2}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}},\end{array}\end{array}$$

where L=L_aL_c. This minimisation problem is regularised for solution boundedness with a parameter λ and given as

$$\begin{array}{}\text{(A8)}& \text{Min}{\mathrm{\Pi }}_E\left(\mathit{\theta }\right)\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{L}\mathit{\theta }-{\mathbf{L}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t{\Vert }_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{\theta }{\Vert }_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

The above minimisation problem is known as Tikhonov regularisation, and λ is referred to as the regularisation parameter. Minimising Eq. (A8) as

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\mathrm{d}{\mathrm{\Pi }}_E}{\mathrm{d}\mathit{\theta }}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\displaystyle \frac{{\mathbf{L}}^{\mathrm{2}}\mathit{\theta }}{\mathrm{2}}}-{\mathbf{LL}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t)+{\mathit{\lambda }}^{\mathrm{2}}\mathit{\theta }=\mathrm{0}\end{array}$$

yields the following quadratic equation in θ:

$$\begin{array}{}\text{(A10)}& \mathit{\theta }=({\displaystyle \frac{{\mathbf{L}}^{\mathrm{2}}}{\mathrm{4}}}+{\mathit{\lambda }}^{\mathrm{2}}\mathbf{I}{)}^{-\mathrm{1}}{\displaystyle \frac{{\mathbf{LL}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t}{\mathrm{2}}},\end{array}$$

where I is the identity matrix of order (N+3).

The choice of regularisation parameter (λ) plays a crucial role in getting an optimal fit for the solution. Based on knowledge about measurement errors, Hansen (2005) proposed two classes for the estimation of λ:

• methods based on knowledge of measurement errors

• methods that do not require details about measurement errors.

In the present scenario, the information regarding the measurement error is unknown; hence class two is used for the current study. In class two, there are three widely used methods (Nieminen et al., 2011): (i) the quasi-optimality criterion, (ii) generalised gross validation (GCV), and (iii) the L-curve method. Compared to the GCV method, the other two methods give a better estimate of λ (Gao et al., 2016). Further, for larger problems, the quasi-optimality method is computationally more expensive than the L-curve method. Owing to this fact, the L-curve method is used here for estimating λ. In the L-curve method, the optimal λ is the one which gives the maximum curvature in the L-curve between the norm of the regularised solution $\mathit{\alpha }\left(\mathit{\lambda }\right)=\Vert {\mathit{\theta }}_{\mathrm{reg}}{\Vert }_{\mathrm{2}}$ and the norm of the residual $\mathit{\beta }\left(\mathit{\lambda }\right)=\Vert \frac{\mathrm{1}}{\mathrm{2}}\mathbf{L}\mathit{\theta }-{\mathbf{L}}_{\mathrm{a}}\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\mathrm{\Delta }t{\Vert }_{\mathrm{2}}$, and the curvature of the L-curve is given by (Nieminen et al., 2011)

$$\begin{array}{}\text{(A11)}& \mathit{\kappa }\left(\mathit{\lambda }\right)={\displaystyle \frac{\ddot{\mathit{\alpha }}\dot{\mathit{\beta }}-\ddot{\mathit{\beta }}\dot{\mathit{\alpha }}}{[{\dot{\mathit{\alpha }}}^{\mathrm{2}}+{\dot{\mathit{\beta }}}^{\mathrm{2}}{]}^{\mathrm{3}/\mathrm{2}}}}.\end{array}$$

Substituting the optimal λ obtained by finding the maximum curvature of Eq. (A11) and $\dot{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\left(t\right)$ obtained from Eq. (3) in Eq. (A10), the regularised torsional displacement (θ) is obtained.

A2 Modal-parameter estimation using the collage method

For a given initial-value problem (IVP),

$$\begin{array}{}\text{(A12)}& x\left(t\right)=f(x,t),x\left(\mathrm{0}\right)=x_{\mathrm{0}},\end{array}$$

that admits a target solution x(t), the associated Picard integral operator T is given by

$$\begin{array}{}\text{(A13)}& \left(\mathrm{T}u\right)\left(t\right)=x_{\mathrm{0}}+\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}f\left(u\right(s),s)\mathrm{d}s.\end{array}$$

The assumptions regarding the parameter estimation problem using the collage method are listed as follows (Deng and Liao, 2009):

1. $x\left(t\right)\in [t_{\mathrm{0}},t_n]$ is a bounded solution, where t₀ and t_n are positive constants satisfying t₀<t_n.

2. $f(u,x,t,{\mathit{\gamma }}_{\mathrm{1}},\mathrm{\cdots },{\mathit{\gamma }}_m)$ is continuous, where ${\mathit{\gamma }}_i,i=\mathrm{1},\mathrm{\cdots },m$ are the unknown modal parameters.

3. The exact solution x(t) of Eq. (A12) uniquely exists.

Here, the unique solution of the considered IVP is given by the fixed point $\overline{u}\left(t\right)$ of this Picard operator (Kunze et al., 2004). Accordingly, the collage distance becomes (x−Tx), and then the optimal solution is the one which minimises the squared ℒ² collage distance (here, ℒ² is 2 norm or square norm of a function). Also, unlike the conventional inverse problem which minimises the approximate error $\mathrm{d}(x-\overline{x})$, the collage method minimises the collage distance d(x,Tx), which is a useful change, as one cannot find $\overline{x}$ for a general T (Kunze et al., 2004). Further, the optimality of the collage distance minimisation is ensured as shown by Kunze et al. (2004). Accordingly, the ℒ² collage distance has the form

$$\begin{array}{}\text{(A14)}& \mathrm{\Delta }=\left(\underset{t_{\mathrm{0}}}{\overset{t_n}{\int }}\right(x\left(t\right)-\left(\mathrm{T}x\right)\left(t\right){)}^{\mathrm{2}}\mathrm{d}t{)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Minimising the squared ℒ² collage distance yields a stationarity condition, $\frac{\mathrm{d}{\mathrm{\Delta }}^{\mathrm{2}}}{\mathrm{d}{\mathit{\gamma }}_{\mathrm{m}}}=\mathrm{0}$, that results in a set of simultaneous linear equations as a function of unknown modal parameters (γ_m). The modal parameters are then obtained by solving that set of linear equations.