Research article 05 Nov 2021
Research article  05 Nov 2021
Identification of wind turbine mainshaft torsional loads from highfrequency SCADA (supervisory control and data acquisition) measurements using an inverseproblem approach
 Department of Wind Energy, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
 Department of Wind Energy, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
Correspondence: W. Dheelibun Remigius (drwp@dtu.dk)
Hide author detailsCorrespondence: W. Dheelibun Remigius (drwp@dtu.dk)
To assess the structural health and remaining useful life of wind turbines within wind farms, the sitespecific structural response and modal parameters of the primary structures are required. In this regard, a novel inverseproblembased methodology is proposed here to identify the dynamic quantities of the drivetrain main shaft, i.e. torsional displacement and coupled stiffness. As a modelbased approach, an inverse problem of a mathematical model concerning the coupledshaft torsional dynamics with highfrequency 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 modelbased 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 sitespecific mainshaft 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 sitespecific mainshaft torsional loads and damageequivalent 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 highfrequency SCADA measurements.
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 sitespecific 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énCosgriff and Kaynia, 2021), and an outputbased 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) timedomainbased methods and (ii) frequencydomainbased methods (Zhang et al., 2010). A recent comprehensive review on various timedomain and frequencydomainbased OMA methods can be found in Zahid et al. (2020). Upon identification of the modal parameters together with the measured structural response, an inverseproblembased technique is employed for the estimation of the sitespecific loads (Pahn et al., 2017). Alternatively, one can use straingaugebased load sensors to measure the sitespecific loads directly (IEC 6140013, 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 SCADAbased (supervisory control and data acquisition) measurements can yield a costeffective 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 inverseproblembased approach is developed here to determine the torsional stiffness and response of the main shaft of a wind turbine, using existing highfrequency SCADA measurements such as the rotor speed and generator speed. This is a costeffective 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 inverseproblem approach is a modelbased 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 highfrequency rotor and generator speed measurements as inputs. Subsequently, the mainshaft torsional displacement is obtained by numerically integrating the shaft torsional velocity. However, numerical integration is based on timemarching algorithms, and the lack of initial conditions makes displacement reconstruction an illposed problem (Hong et al., 2008). Since it is an illposed problem, influence of measurement noise will be amplified during the timemarching procedure, resulting in an erroneous displacement. This inaccurate displacement also leads to drastic errors in the system modalparameter estimation. Hence, one needs to go for regularisation techniques to smoothen the reconstructed displacement. Hansen (2005) discussed the nature of various illposed 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 infiniteimpulseresponse (IIR) and finiteimpulseresponse (FIR) filters have several drawbacks when reconstructing the lowfrequency 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 lowfrequency dominant structures. Hence, the same has been employed in the present work. Tikhonov regularisation minimises the error using the leastsquare 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 modelbased system identification technique – is used (Kunze and Vrscay, 1999; Groetsch, 1993). This method has been successfully applied for the system identification in various differentialequationbased 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 inverseproblembased methods (Deng and Liao, 2009). Further, the modelbased collage method is simple, easy to implement, and computationally inexpensive as compared to the outputbased OMA methods.
The estimated shaft torsional stiffness and displacement are further used to identify the sitespecific 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 sitespecific load can be determined without requiring additional sensors or information from the wind turbine manufacturer. The mainshaft torsional load affects the fatigue performance of other drivetrain components such as the gearbox and planetary bearings (Dong et al., 2012; GallegoCalderon and Natarajan, 2015; Ding et al., 2018). Hence, the same sitespecific 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.
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 inverseproblem algorithm and the estimated shaft torsional displacement θ, and torsional stiffness K will be utilised for the shaft torsionalload estimation. For this purpose, a twomass model (refer to Fig. 1) (Boukhezzar et al., 2007; Girsang et al., 2013; Berglind et al., 2015) that governs the mainshaft 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 lowspeed 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 twomass model (Girsang et al., 2013).
Here, J_{r} represents the inertia of the rotor; J_{g} represents the collective inertias of the highspeed 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 forwardproblem 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 inverseproblem 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.
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 numericalintegration 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.
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
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{\u0303}}{\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 numericalintegration 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 numericalintegration result.
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 modelbased 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 ℒ^{2} collage distance (refer to Appendix A2) for Eq. (2) becomes
where ${\mathit{\omega}}_{{\mathrm{g}}_{\mathrm{0}}}$ is the generator speed at time, t=t_{0}. 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 mainshaft 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) VestasV52 (Vestas, 2020). For facilitating a comparison of the mainshaft 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 inhouse 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 hydroeffects (in the case of offshore turbines) on the structure. Control of the turbine is performed through the dynamic link libraries (DLLs). Using the mainshaft torsionalload 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 torsionalstiffness 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.
Aeroelastic simulations are performed for the Vestas V52 turbine corresponding to the design load case (DLC 1.2) (IEC 614001, 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.
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
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 lowspeed side is as follows:
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
Finally, the static and dynamic components are superimposed to get the actual displacement (θ). This actual displacement will then be used to estimate the sitespecific mainshaft torsional load as shown by Eq. (9).
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 timeseries variation and the dominant frequency dynamics (lowfrequency 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 highperformance computing cluster of DTU. The cluster has 320 nodes in total, with each node consisting of two Intel Xeon E52680 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.
Upon estimating the torsional load, the torsional damageequivalent load (DEL) at each mean wind speed is calculated using the following equation:
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 zeromean 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 mainshaft 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 mainshaft 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.
Upon verifying the proposed method, the drivetrain mainshaft 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 lowspeed 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 wakefree condition between 280 and 320^{∘} wind directions.
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 Lcurve 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 sitespecific 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 mainshaft 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 torsionalload time series, which may be used an input to gearbox design tools to predict the loading within the gearbox.
A novel inverseproblembased approach is developed for estimating the mainshaft torsional displacement and stiffness by using highfrequency SCADA measurements. A mathematical model describing the coupledshaft 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 mainshaft coupled torsional stiffness. The estimated mainshaft quantities are then used to identify the mainshaft sitespecific torsional load. The proposed formulation is successfully verified for the mainshaft torsional loads with the aeroelastic simulation of the Vestas V52 turbine. Upon verification, the methodology is extended to identify the sitespecific mainshaft 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 torsionalload identification. Since the sitespecific 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 highfrequency SCADA measurements for the estimation of the mainshaft torsional load and DEL.
A1 Displacement reconstruction using Tikhonov regularisation
By definition, the velocity $\dot{\mathit{\theta}}$ is expressed as
where $\dot{\stackrel{\mathrm{\u0303}}{\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 leastsquare sense, the following minimisation problem has to be solved:
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):
where Δt is the time interval of the discretisation and L_{a} is the diagonal weighing matrix of order (N+1) as
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 centraldifference rule and written in matrix form as
where the centraldifference matrix L_{c} of size $(N+\mathrm{1})\times (N+\mathrm{3})$ and the displacement vector θ of size (N+3) are given as
In the finitedifference 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 finitedifference method (Lapidus and Pinder, 2011). Substitution of Eq. (A5) into Eq. (A3) leads to
where L=L_{a}L_{c}. This minimisation problem is regularised for solution boundedness with a parameter λ and given as
The above minimisation problem is known as Tikhonov regularisation, and λ is referred to as the regularisation parameter. Minimising Eq. (A8) as
yields the following quadratic equation in θ:
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 quasioptimality criterion, (ii) generalised gross validation (GCV), and (iii) the Lcurve method. Compared to the GCV method, the other two methods give a better estimate of λ (Gao et al., 2016). Further, for larger problems, the quasioptimality method is computationally more expensive than the Lcurve method. Owing to this fact, the Lcurve method is used here for estimating λ. In the Lcurve method, the optimal λ is the one which gives the maximum curvature in the Lcurve 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{\u0303}}{\mathit{\theta}}}\mathrm{\Delta}t{\Vert}_{\mathrm{2}}$, and the curvature of the Lcurve is given by (Nieminen et al., 2011)
Substituting the optimal λ obtained by finding the maximum curvature of Eq. (A11) and $\dot{\stackrel{\mathrm{\u0303}}{\mathit{\theta}}}\left(t\right)$ obtained from Eq. (3) in Eq. (A10), the regularised torsional displacement (θ) is obtained.
A2 Modalparameter estimation using the collage method
For a given initialvalue problem (IVP),
that admits a target solution x(t), the associated Picard integral operator T is given by
The assumptions regarding the parameter estimation problem using the collage method are listed as follows (Deng and Liao, 2009):

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

$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.

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 ℒ^{2} collage distance (here, ℒ^{2} 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 ℒ^{2} collage distance has the form
Minimising the squared ℒ^{2} 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.
The code regarding the mathematical models developed in the article can be accessed at https://github.com/dheelibun/EstimationofmainshafttorsionalmodalparametersfromhighfrequencySCADA#readme (last access: 4 November 2021) and https://doi.org/10.5281/zenodo.5643220 (Remigius, 2021).
The opensource DTU 10 MW wind turbine aeroelastic model can be accessed at https://rwt.windenergy.dtu.dk/dtu10mw/dtu10mwrwt (last access: 4 November 2021; Zahle, 2018).
Vestas V52 wind turbine SCADA data and its other parameters are not publicly available. However, Vestas v52 SCADA data can be requested by signing a nondisclosure agreement.
WDR conceived the methodology; completed the formal analysis, investigation, and validation of the project; and wrote the original draft of the paper. AN conceived the original idea, developed the scientific methods, reviewed and edited the paper, and supervised the project.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work has been funded by the Energy Technology Development and Demonstration Program (EUDP) of the Energistyrelsen, Denmark (LIFEWIND project; grant no. 6401705114). Investigations were carried out at the Department of Wind Energy, Technical University of Denmark.
This paper was edited by Amir R. Nejad and reviewed by Edward Hart and one anonymous referee.
Aarden, P.: System Identification of the 2B6 Wind Turbine: A regularised PBSIDopt approach, Master's thesis, Delft University of Technology, Delft, the Netherlands, 2017. a
Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L. C., Natarajan, A., and Hansen, M.: Description of the DTU 10 MW reference wind turbine, Tech. Rep. I0092, DTU Wind Energy, Roskilde, Denmark, 2013. a
Bangji, Z., Shouyu, Z., Qingxi, X., and Nong, Z.: Load identification of virtual iteration based on Tikhonov regularization and model reduction, Hong Kong Journal of Social Sciences, 44, 53–59, 2017. a
Berglind, J. J. B., Wisniewski, R., and Soltani, M.: Fatigue load modeling and control for wind turbines based on hysteresis operators, in: 2015 American Control Conference (ACC), IEEE, Chicago, IL, USA, 3721–3727, 2015. a
Boukhezzar, B., Lupu, L., Siguerdidjane, H., and Hand, M.: Multivariable control strategy for variable speed, variable pitch wind turbines, Renew. Energy, 32, 1273–1287, 2007. a, b
Brandt, A. and Brincker, R.: Integrating time signals in frequency domain–Comparison with time domain integration, Measurement, 58, 511–519, 2014. a
Deng, X. and Liao, Q.: Parameter estimation for partial differential equations by collagebased numerical approximation, Mathematical Problems in Engineering, 2009, 510934, https://doi.org/10.1155/2009/510934, 2009. a, b
Deng, X., Wang, B., and Long, G.: The Picard contraction mapping method for the parameter inversion of reactiondiffusion systems, Comput. Math. Appl., 56, 2347–2355, 2008. a
Ding, F., Tian, Z., Zhao, F., and Xu, H.: An integrated approach for wind turbine gearbox fatigue life prediction considering instantaneously varying load conditions, Renew. Energy, 129, 260–270, 2018. a
Dong, W., Xing, Y., and Moan, T.: Time domain modeling and analysis of dynamic gear contact force in a wind turbine gearbox with respect to fatigue assessment, Energies, 5, 4350–4371, 2012. a
GallegoCalderon, J. and Natarajan, A.: Assessment of wind turbine drivetrain fatigue loads under torsional excitation, Eng. Struct., 103, 189–202, 2015. a
Gao, W., Yu, K., and Wu, Y.: A new method for optimal regularization parameter determination in the inverse problem of load identification, Shock and Vibration, 7328 969–1–16, 2016. a
Girsang, I. P., Dhupia, J. S., Muljadi, E., Singh, M., and Pao, L. Y.: Gearbox and drivetrain models to study dynamic effects of modern wind turbines, in: 2013 IEEE Energy Conversion Congress and Exposition, 874–881, https://doi.org/10.1109/ECCE.2013.6646795, 2013. a, b
Groetsch, C. W.: Inverse problems in the mathematical sciences, vol. 52, Springer, Wiesbaden, Germany, 1993. a, b
Hansen, P. C.: Rankdeficient and discrete illposed problems: Numerical aspects of linear inversion, vol. 4, SIAM, Philadelphia, USA, 2005. a, b, c
Hào, D. N. and Quyen, T. N. T.: Convergence rates for Tikhonov regularization of a twocoefficient identification problem in an elliptic boundary value problem, Numerische Mathematik, 120, 45–77, 2012. a
Hedges, S. B. and Shah, P.: Comparison of mode estimation methods and application in molecular clock analysis, BMC Bioinformatics, 4, 31–1–11, 2003. a
Hong, Y. H., Park, H. W., and Lee, H. S.: A regularization scheme for displacement reconstruction using acceleration data measured from structures, in: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 2008, International Society for Optics and Photonics, 6932, 693228, 2008. a, b, c
IEC 614001: Wind energy generation systems – Part 1: Design requirements, Standard, International Electrotechnical Commission, Geneva, Switzerland, 2019. a
IEC 6140013: Wind turbines – Part 13: Measurement of mechanical loads, Standard, International Electrotechnical Commission, Geneva, Switzerland, 2016. a
Koukoura, C., Natarajan, A., and Vesth, A.: Identification of support structure damping of a full scale offshore wind turbine in normal operation, Renewable Energy, 81, 882–895, 2015. a
Kunze, H. and La Torre, D.: Computational Aspects of Solving Inverse Problems for Elliptic PDEs on Perforated Domains Using the Collage Method, in: Mathematical and Computational Approaches in Advancing Modern Science and Engineering, Springer, Cham, Germany, 113–120, 2016. a
Kunze, H., La Torre, D., and Vrscay, E.: A generalized collage method based upon the Lax–Milgram functional for solving boundary value inverse problems, Nonlinear Analysis: Theory, Methods & Applications, 71, e1337–e1343, 2009. a, b
Kunze, H. E. and Vrscay, E. R.: Solving inverse problems for ordinary differential equations using the Picard contraction mapping, Inverse Problems, 15, 745–770, 1999. a, b, c
Kunze, H. E., Hicken, J. E., and Vrscay, E. R.: Inverse problems for ODEs using contraction maps and suboptimality of the “collage method”, Inverse Problems, 20, 977–991, 2004. a, b, c
Lapidus, L. and Pinder, G. F.: Numerical solution of partial differential equations in science and engineering, John Wiley & Sons, New York, USA, 2011. a
Larsen, T. J. and Hansen, A. M.: How 2 HAWC2, the user's manual, Tech. rep., Risø National Laboratory, Technical University of Denmark, Roskilde, Denmark, 2007. a
Lee, H. S., Hong, Y. H., and Park, H. W.: Design of an FIR filter for the displacement reconstruction using measured acceleration in lowfrequency dominant structures, Int. J. Numer. Meth. Eng., 82, 403–434, 2010. a, b, c, d, e
Nieminen, T., Kangas, J., and Kettunen, L.: Use of Tikhonov regularization to improve the accuracy of position estimates in inertial navigation, International Journal of Navigation and Observation, 2011, 450269, https://doi.org/10.1155/2011/450269, 2011. a, b, c
NorénCosgriff, K. and Kaynia, A. M.: Estimation of natural frequencies and damping using dynamic field data from an offshore wind turbine, Marine Struct., 76, 102915, 2021. a
Pagitsch, M., Jacobs, G., and Bosse, D.: Remaining Useful Life Determination for Wind Turbines, J. Phys.Conf. Ser., 1452, 012052, https://doi.org/10.1088/17426596/1452/1/012052, 2020. a
Pahn, T.: Inverse load calculation for offshore wind turbines, PhD thesis, Gottfried Wilhelm Leibniz Universität Hannover, Hannover, 2013. a
Pahn, T., Rolfes, R., and Jonkman, J.: Inverse load calculation procedure for offshore wind turbines and application to a 5MW wind turbine support structure, Wind Energy, 20, 1171–1186, 2017. a, b, c
Qihe, L.: Integration of vibration acceleration signal based on labview, J. Phys.Conf. Ser., 1345, 042067, https://doi.org/10.1088/17426596/1345/4/042067, 2019. a
Remigius, W. D.: Estimation of main shaft torsional modal parameters from highfrequency SCADA data, in: Identification of wind turbine mainshaft torsional loads from highfrequency SCADA (supervisory control and data acquisition) measurements using an inverseproblem approach, v1.0, 6, 1–12, Zenodo [code], https://doi.org/10.5281/zenodo.5643220, 2021. a
Ronasi, H., Johansson, H., and Larsson, F.: A numerical framework for load identification and regularization with application to rolling disc problem, Comput. Struct., 89, 38–47, 2011. a
Tikhonov, A. N.: Solution of incorrectly formulated problems and the regularization method, Soviet Mathematics – Doklady, 4, 1035–1038, 1963. a, b
Veldkamp, H. F.: Chances in wind energy: A probabilistic approach to wind turbine fatigue design, PhD thesis, Faculty of Mechanical Maritime and Materials Engineering, TU Delft, 2006. a, b, c, d
Vestas: Vestas V52, available at: https://en.windturbinemodels.com/turbines/71vestasv52, last access: 8 April 2020. a
Wang, T., Celik, O., Catbas, F. N., and Zhang, L. M.: A frequency and spatial domain decomposition method for operational strain modal analysis and its application, Eng. Struct., 114, 104–112, 2016. a
Zahid, F. B., Ong, Z. C., and Khoo, S. Y.: A review of operational modal analysis techniques for inservice modal identification, J. Braz. Soc. Mech. Sci., 42, 1–18, 2020. a
Zahle, F.: dtu10mwrwt, DTU [code], available at: https://rwt.windenergy.dtu.dk/dtu10mw/dtu10mwrwt (last access: 4 November 2021), 2018. a
Zhang, L., Wang, T., and Tamura, Y.: A frequency–spatial domain decomposition (FSDD) method for operational modal analysis, Mechanical systems and signal processing, 24, 1227–1239, 2010. a
Ziegler, L., Gonzalez, E., Rubert, T., Smolka, U., and Melero, J. J.: Lifetime extension of onshore wind turbines: A review covering Germany, Spain, Denmark, and the UK, Renewable and Sustainable Energy Reviews, 82, 1261–1271, 2018. a, b