Digital twin modeling for predictive maintenance of gearboxes in floating offshore wind turbine drivetrains

This paper presents a multi-degree of freedom torsional model of drivetrain system as the digital twin model for monitoring the remaining useful lifetime of the drivetrain components. An algorithm is proposed for the model identification, which receives the torsional response and estimated values of rotor and generator torques, and calculates the drivetrain dynamic properties, e.g. eigenvalues, and torsional model parameters. The applications of this model in prediction of gearbox remaining useful lifetime is discussed. The proposed method is computationally fast, and can be implemented by integrating with the current turbine control and monitoring system without a need for a new system and sensors installation. A test case, using 5 MW reference drivetrain, has been demonstrated.


Introduction
Power train system including rotor, main bearings, gearbox, generator and power converter accounts for 57% of turbine total failures and 65% of turbine total downtime [1]. The consequent costs are expected to be higher in floating wind turbines (FWT) due to costly and slow-going marine operations. In FWTs, the power train system is exposed to a wider range of excitations than bottom-fixed turbines due to synergistic impacts of wind, wake, currents and structural Farid K. Moghadam farid.k.moghadam@ntnu.no 1 Norwegian University of Science and Technology, Trondheim, Norway motions-induced vibrations [2]. Predictive maintenance is proposed in the literature as a support tool for the time base maintenance to reduce the unexpected shutdowns and the consequent downtimes while it can also help to optimize the inspection intervals. The motivation of this research is increasing the wind turbine (WT) availability by performing predictive maintenance of the drivetrain system through the digital twin models. Digital twins are highly accurate but computationally fast models of the system which can be run in-line with the real machine. The model should be able to update itself by the online operational measurements, to capture the physical variations in the system over the time [3]. In this paper, digital twin is proposed as a tool for lifetime monitoring of the drivetrain system, and more specifically the gearbox, by using the torsional vibration measurements, which K can be included in the farm Supervisory Control and Data Acquisition (SCADA) existing system by using higher resolution encoders. The input data may also be provided by additional installations e.g. proximity sensors and angular accelerometers. Correa et al. [4] performs a review on the drivetrain components condition monitoring techniques and faults prediction in general and also by using the SCADA data. Digital twin is a proven technology which is used by Siemens Gamesa for prediction of drivetrain loads and subsequently improving the drivetrain design [5]. The application of digital twin for monitoring of the drivetrain lifetime may be looked as an expensive solution. However, in this work, it is shown that by using innovative approaches based on analysis of torsional measurements it is possible to offer computationally-fast solutions for monitoring the lifetime of the drivetrain critical components with a potential to be integrated with the currently available control and monitoring system with minimum additional cost. Digital twin in this context means the combination of model, online measurements and remaining useful lifetime (RUL) model as defined and suggested by [6].
Rebouças et al. [7] suggests the 14-DOF torsional model as a fair compromise between complexity and accuracy when the structural integrity of gearbox is taken into consideration. The structural integrity is mainly measured in terms of contact fatigue stresses of the gear pairs, so that this model can be used as the indicator for the gears contact loads. For this purpose, the applicability of a 14-DOF lumped-parameter torsional model of the drivetrain system as the digital twin model for monitoring the degradation in a high-speed drivetrain technology with a three-stages gearbox which consists of two planetary and one parallel stages is investigated. The algorithm for estimation of digital twin model parameters is based on a data-driven approach which operates based on the real-time drivetrain torsional measurements, and the estimated dynamic properties from the torsional response. The algorithm employed for estimation of dynamic properties, namely natural frequencies and damping from the torsional response is according to the modal estimation approach proposed by [8]. In order to estimate the moment of inertia of the components in the equivalent model, an optimization problem is defined based on the least square error between the digital twin model and the online measurements consisting of torsional response and estimated input torque. After all, the stiffness parameters of the model are estimated by using the equations of natural frequencies as a function of equivalent model parameters (stiffness and inertia). It is worth noting that the equivalent lumped parameter model of drivetrain in the proposed digital twin approach is linear. However, the eigenfrequencies of the drivetrain are nonlinear functions of equivalent model parameters. Since the assumption of accessibility of all the system dynamic states is opti-mistic, the problem is extend to a more general class of problem with a partial knowledge about the system modes. The parameters of the digital twin model are updated by using the real-time operational data. The most prevalent failures of the large gears in WT drivetrain systems are due to gear tooth root bending and pitting fatigue damage [9]. For monitoring the RUL of the gears in this work, the degradation of the gears due to pitting fatigue damage and the gear pairs contact loads are taken into consideration. The contact loads on the gear pairs are estimated by using the load observers designed for the gear pairs based on the real-time estimated digital twin model and torsional measurements. The employed model-based degradation model is based on stress-life method. The online estimated contact stress feeds the time-domain cycle counting approach based on rainflow method for online estimation of fatigue cycles, and then for estimation of fatigue damage for each gear pair by using the Miner's rule. For all the simulation cases, the floating 5 MW NREL WT with a spar support substructure as a common way of realization of FWT is used. The reason for the latter is the higher motivation in FWT for predictive maintenance.
This paper is the proof of concept for the idea of near real-time estimation of the load and subsequently the residual life in the drivetrain components by using computationally efficient equivalent models estimated from torsional measurements supported by real-time values of torsional response. The proposed algorithm is able to estimate the equivalent model parameters by using only few measurement samples which feed a robust linear regression estimator designed to estimate the equivalent model parameters. On this basis, the main contributions of this work are: Proposing the 14-DOF equivalent lumped-parameter model as the digital twin of the drivetrain system for remaining useful lifetime monitoring of the gearbox, Proposing an algorithm for the near real-time estimation of the equivalent model parameters by using the torsional measurements, Designing computationally inexpensive load observers for the near real-time estimation of the contact load and stress on the gears in planetary and parallel stages of the drivetrain gearbox, by using equivalent model parameters and real-time torsional responses, Proposing a physics-based degradation model for near real-time estimation of the residual life of the gears in the drivetrain gearbox by using the online estimated equivalent model, real-time measurements, designed load observers and contact stress estimation approach, Simulation studies to evaluate the proposed digital twin approach for estimation of fatigue damage of the gears, and validation of the results by high-fidelity simulation models.
K Fig. 1 Representation of the torsional model [7] 2 Methodology

Test case
NREL 5 MW reference wind turbine [10] with a spar floating support substructure is selected for this study. The wind turbine specification and the overall characteristics of the floating platform is obtained from [11]. This model is able to capture the global dynamics of spar floating wind turbine from the interaction with the environmental loads. The 5 MW reference drivetrain [12] is employed in this study.

Global simulation and estimation of drivetrain loads
The decoupled simulation approach is used for the drivetrain studies in this work. The latter means that the turbine global simulation results, namely the rotor aerodynamic load and the bedplate motions are applied to the detailed local drivetrain model in a secondary software. Then the drivetrain local load effects are obtained for further postprocessing aimed at monitoring the drivetrain remaining lifetime. More details on the utilized decoupled approach and the global simulation model is available in [11]. The wind is turbulent based on Kaimal distribution. The turbulence intensity at hub height I (-) is assumed to be 0.14, for all the wind speeds, according to IEC 61400-1 class B turbines. The wave is modelled stochastic by two parameters, namely significant wave height H s .m/and peak period T p .s/ in global analysis.

Drivetrain torsional model as the digital twin model for predictive maintenance
The NREL 5 MW reference drivetrain proposed in [12] is considered in this study, considering only its torsional behavior. Fig. 1 contains a schematic representation of the model used in this work. The gearbox is assumed to be clamped to the nacelle's bedplate by rigid torque arm supports, leading to no rotation of the ring gears. This model is based on [13], where the resonant behavior of a single-stage planetary gearbox is discussed. The equations of motion are shown in Eq. (1).
K In these equations, Â k and T k (k 2 fR; 1; 2; 3; Gg) represent torsional coordinate and external torques, respectively. The indexes R and G are associated with rotor and generator variables, respectively. The sub-vectors Â i (i 2 f1; 2; 3g) refer to the gear stages, see Eqs. (3) and (6). The global inertia and stiffness matrices J and K are obtained from their local counterparts. The damping is defined using Rayleigh parameters, being proportional to the inertia and stiffness matrices J and K, i.e D =˛J +ˇK.
The local matrices for the planetary stages can be seen in Eqs. (3) to (5) and are based on the matrices used in [13], while the local matrices for the parallel stage, shown in Eq. (6), were obtained by isolating the sun-planet gear pair of the planetary stage. In these equations, Â, m, k, and r, represent the angular displacement, mass, stiffness and base radius respectively. The indices c; p; s; P , and W represent the carrier, planet and sun elements of the i th stage, respectively. The mesh stiffness between sun-planet, ring-planet and pinion-wheel gears is represented by k sp , k rp , and k P W , respectively. In this work, the mesh stiffness between different gear pairs is calculated according to ISO 6336-2 [14], which takes a constant value. The center distance for planetary stages is represented by a w [7]. Shafts are used to connect the adjacent inertia elements, such as the rotor, generator, and gear stages. The shafts were modelled using finite element theory leading to the shaft's inertia and stiffness matrices J S and K S in Eq. (7). Assembling the global matrices from the local matrices is illustrated in Fig. 2. In Eq.
The generator torque is realized by a proportional-integral (PI) speed controller, defined as where the error e is given as the difference between reference and measured generator speeds ! G and P Â G , respectively. K P = 2200 is the proportional gain and K I = 220 is the integral gain.

Estimation of model parameters from torsional response
Drivetrain modal estimation The torsional response residual function between the inertia j m and j n from the point m is defined as em ;n ,˛m − u m;n˛n ; for m; n 2 f1; :::; 14g ; with˛is the angular acceleration, and u m;n is the relative gear ratio between j m and j n to make them in the same coordinate. Gear ratio u m;n as per definition is N m N n , where N m and N n are the speeds at m th and n th inertias. The theoretical explanation of estimation of natural frequencies from the frequency spectrum of em ;n is provided by [8]. The different residual functions can be defined by using the response obtained from the different locations of drivetrain, which present different levels of visibility of the drivetrain torsional modes.

Estimation of moment of inertia matrix
By summing the moments on each of the inertias of the lumped parameter model yields 14 equations of the form for i = .1; :::; 14/ : In the matrix form, these set of equations can be written as K Forsch Ingenieurwes Assuming that the load and response time series are known, the parameter estimation turns to the minimization of the L2 norm of error. The error function is defined by The least square estimator is defined by where J is a diagonal matrix. K and D are nondiagonal but symmetric matrices. The matrices D and K are not full rank. The latter causes some computational difficulty for the above quadratic matrix optimization problem, so that there is no guarantee for the convergence. In order to remove the coupling between these equations due to the K and D terms, and to reduce the computational complexity by reduction of the number of variables, the equivalent scalar optimization problem is constructed by the sum of the individual equations of each inertia. The latter leads to the following error function in terms of the variables j i with rotor as the reference of the rotary coordinate.
e .t/ = j 1˛1 .t/ + ::: + a 1;i j i˛i .t/ + :::+ where˛i is the time series of angular acceleration at the i th body. Rot and Gen are the time series of the rotor and generator torques respectively. The sign of a 1;i is determined based on the direction of rotation of j i . The latter leads to the following quadratic scalar optimization problem as This estimator is robust to the measurement noises, and can provide a good approximation even with less than 14 input data samples (underdetermined case). For the case of more than 14 samples (overdetermined case), this estimator helps to obtain more accurate estimation than solving the linear equations, when the input measurements are subject to independent identically distributed (i.i. d.) Gaussian noise. In other words, the total LS technique is able to correct the system with minimum perturbation [15]. The above convex optimization problem is numerically solved by Matlab CVX and the global optimizer j LS = fj 1 ; :::; j 14 g is estimated.
Estimation of stiffness matrix The undamped torsional frequencies of the system are the nonlinear function of inertia and stiffness as ! i .for i = 1; ::: By using the estimated natural frequencies obtained from the modal estimation approach together with the estimated inertia matrix J from the LS optimization problem, the stiffness matrix K is the root of g i which is defined by the following nonlinear equation as .for i = 1; :::; 14/: In general, there is not a unique matrix K from the above equation for the known set of eigenvalues i = f! 1 2 ; :::; ! 14 2 g of the matrix −J −1 K. However, by imposing the sparsity and symmetricity to matrix K from the lumped model, it is possible to calculate the unique matrix K numerically by using Matlab fsolve solver. The latter also helps to reduce commutation cost of this matrix algebraic equation by reducing the number of variables from N 2 to N. The matrix J −1 K is not symmetric in the general case which may give the sense that there are multiple answers for K from this equation. However the fact that −J −1 K always has positive eigenvalues (it is positive definite), brings us to the believe that this matrix is a small perturbation of a symmetric matrix with positive eigenvalues. Small perturbation keeps the eigenvalues positive [16].
The usual condition for the estimation problem is more restrictive. In other words, it is possible that only some of the eigenfrequencies of the drivetrain system can be estimated by employing the aforedescribed modal estimation approach, especially the higher eigenfrequencies which are excited with a lower energy of the input torque. In this case, the matrix K can still be estimated by using the following optimization problem in terms of the first n eigenfrequencies as defined by the following least square error estimator: with ƒ is the variable of this problem which is a function of the unknown variable K as ƒ = −J −1 K. Also b k LS is the set of nonzero elements of matrix K which are estimated by the above nonlinear matrix optimization problem. The sign of the elements of k are forced in the optimization problem. n is the set of n .n 2 f1; :::; 14g/ smallest magnitude eigenvalues which are known from the modal estimation, n = f! 1 2 ; :::; ! n 2 g. eig.ƒ; n/ is the set of n ..n 2 f1; :::; 14g/) smallest magnitude eigenvalues de- where S is the sparsity of matrix ƒ. The positive definiteness and sparsity of ƒ are the nonlinear constraints which are imposed to this problem. For the set of positive semidefinite matrices, this problem is convex and the solution is the global optimizer. However, ƒ is not symmetric in general so that the definition of the problem is nonconvex for the numerical solvers and convex optimization tools are not able to numerically solve the problem. For this purpose, Matlab fmincon solver as a powerful tool for the general class of nonlinear nonconvex problems is used.

Estimation of gearbox loads, and defining degradation model
Estimating contact stresses and loads at each stage The free-body diagram for the gears in planetary and parallel stages can be seen in Fig. 3. For wind turbine gearboxes, the input and output torques for planetary (parallel) stages are the carrier (pinion) torque T C .T P / and the sun (wheel) torque T s .T W /, respectively. From Fig. 3, one obtains Eqs. (20) and (21) as Planetary stage: Parallel stage: One can obtain the relationships between the input and output torques by eliminating the internal forces between elements. Speed relations are time differentiated which enables writing the contribution of inertial torques by using a single torsional acceleration and equivalent mass moment of inertia. For the planetary stage, one obtains with J EQ as the equivalent mass moment of inertia, having contributions from the sun, N p planet gears and planet carrier. Similarly, for the parallel stage one has with J EQ having contributions from both pinion and wheel gears. The term multiplying the torque at the pinion is the gear ratio for parallel gears. Torque transfer between stages is made via shafts, and can be estimated by K where J S R Â S is the inertial torque from the shaft, [17]. Therefore, for the NREL 5 MW drivetrain, one should have The input torques estimated above can be used to estimate the stresses at the different gear stages as described below.
Gear contact stresses are analyzed in this work following ISO 6336-2:2019 [14]. According to this standard, the contact stresses are defined as where u is the gear ratio of the pair, d 1 and b are reference diameter and face width of the pinion, respectively, T i is the input torque. The other parameters account for different aspects of the problem, such as contact relations Z BD and ratios Z , material properties Z E , helix angle Zˇ, mesh load K , gear speed K v , load distributions K Hˇi and K H˛i . These factors are discussed in [15] and references therein. Eq. (26) can be rewritten in a compact form as where C represents the K and Z design parameters mentioned above. This parameter, unknown to the maintenance engineer, can be roughly estimated from nominal conditions as C = HN = p T iN . Therefore, Eq. (27) turns to where H i is defined in terms of nominal contact stress HN and torque T iN , which have a clear physical meaning. Additionally, one can also expand Eq. (28) using Taylor series around the nominal torque T iN as This formula gives a polynomial relationship between stress and load deviations. The accuracy of this expression A limitation of this derivation is that C depends on parameters that may change with T i , but these variations should be negligible.

Remaining useful lifetime estimation
The real-time accumulated damage is estimated as follows. First, the timevarying gear transmitted loads of the gear pairs of three gearbox stages are estimated by using the torque observers designed as explained in the previous part. The latter is based on the real-time aerodynamic and generator torques, and the time-varying digital twin model parameters. Then the gear tooth surface pitting stress is estimated for the different gears as explained on above. The number of gear tooth contact stress cycles at different stress levels is counted by using the time-domain rainflow cycle counting approach [18]. The outputs are the amplitude stress level s , and the number of stress cycles at s for s = .1; :::; S /. To consider the influence of nonzero mean stress level, Goodman rule is employed to calculate the effective stress (the equivalent zero mean alternating stress) by [19], where m and u are the mean stress and material yield strength, respectively. The accumulated pitting damage for the data block t with S different stress levels e s .s 2 f1; :::; S g) is calculated by using Miner's rule as where n s is the number of cycles at the stress level e s and N s is the number of cycles to yield at stress level e s , where N s = k. e s / −m . The total absolute online accumulated damage will then be calculated by where T stands for the last data block which represents the current time. This method can also be used for estimation of relative damage between different operational periods over K the time, to give an insight on variations in degradation between different operational periods.

Aerodynamic and generator torques, and drivetrain torsional response
The input aerodynamic torque obtained from SIMO-RIFLEX-AeroDyn and the performance of the designed controller in controlling the generator torque under variable input torque to set the speed on the shaft is shown in Fig. 4. The drivetrain model responses are the angular displacement, velocity and acceleration in the different bodies in the described 14-DOF model. As an example, the angular velocity responses of the rotor and generator are shown in Fig. 4. For demonstration purposes the torques a b and response are scaled with the gear ratio. The turbine operation is assumed to be near the rated operation.

Estimation of dynamic properties and digital twin model parameters
The undamped frequency modes of the 14-DOF model for the healthy system are listed in the Table 1. The torsional response error function can be defined between different bodies in the drivetrain model. For example, the undamped natural frequencies estimated from the angular acceleration error function for the 2nd gearbox stage is shown in Fig. 5. The main feature of angular acceleration compared to angular velocity and displacement is the amplification of the higher natural frequencies in the response. By defining the angular acceleration error functions in terms of different pairs of bodies in the model, it is possible to estimate the different drivetrain modes. However, the significance  0  2  120  153  208  567  716  716  971  1224  1232  1232  1384  1943  Fault sun   stage1   0  2  120  153  208  565  700  700  929  1224  1232  1232  1384  1943  Fault sun   stage2   0  2  119  153  208  567  716  716  971  1220  1220  1224  1383  1850  Fault   ring  stage1   0  2  120  153  208  541  696  696  968  1224  1232  1232  1384  1943  Fault   ring  stage2   0  2  120  153  208  567  716  716  971  1181  1181  1224  1314  1941 of modes is different in the error functions defined between different bodies. As it can be seen in Fig. 5, the 6th and 7th modes are not visible in the error function of the 2nd stage, but they can be observed in the error function defined between rotor and generator bodies. The estimated natural frequencies by using the employed modal estimation approach based on angular acceleration error function are listed in the Table 1. The maximum relative error of estimation is less than 1.5%. In practice, the torsional response in commercially available turbines is only accessible by the encoders placed on rotor and possibly generator shafts, so that some of the modes may not be observable by only using the error function between rotor and generator. The influence of pitting on parameters of drivetrain equivalent model is usually represented by a decrease on the mesh stiffness of the faulty gear pair [20,21]. The NREL 5 MW drivetrain presented in Sect. 2.3 is used to illustrate this behavior, assuming that the early occurrence of pitting at the sun gear can lead to a 10% decrease in the sun-planet mesh stiffness. The 10% reduction of the faulty gear mesh stiffness assumed for modelling a very early stage fault might be high, but this assumption is made to demonstrate more clearly the influence of the fault on the variations of drivetrain dynamic properties. The results of the mentioned fault in the first and second gear stages can be seen in the second and third rows of Table 2, where one can see that pitting in the first stage sun gear affected mainly the 7th, 8th and 9th resonances, whereas pitting in the second stage sun gear affected mainly the 14th resonances. The results of fault simulation reported in Table 2 show that a 10% reduction in the sun-planet mesh stiffness of the first gear stage, results in around 2.5%, 2.5% and 4.5% reduction of 7th, 8th and 9th drivetrain natural frequencies, respectively, and a 10% reduction in the sun-planet mesh stiffness of the second gear stage, results in about 5% reduction of the highest natural frequency of drivetrain. The influence of pitting fault in the ring of first and second gear stages on the drivetrain natural frequencies are listed respectively in the fourth and fifth rows of Table 2. A 10% reduction in the ring-planet mesh stiffness of the first gear stage, results in around 5%, 3% and 3% reduction of 6th, 7th and 8th natural frequencies, and a 10% reduction in the ring-planet mesh stiffness of the second gear stage, results in the 3.5%, 4% and 5% reduction of 10th, 11th and 13th natural frequencies, respectively. Since the relative error of estimating the natural frequencies can be as high as 1.5%, variations of the natural frequencies which are less than 1.5% cannot be used as the indicator of fault. More detailed models of pitting fault can be engaged to capture more precisely the dynamics of fault [22] and consequent influence on the drivetrain dynamic properties, which is not the scope of this work. The influence of drivetrain faults at system-level on dynamic properties is analytically studied by [23]. Monitoring the variations in the drivetrain system natural frequencies can support the fault detection of the drivetrain at component-level, e.g. faults in the gears of the gearbox gear stages, which is also not the main scope of this work.
The actual values of moment of inertia of the different bodies in the 14-DOF model and the performance of the designed LS estimator in estimation of those parameters from the torsional measurements based on the theory elaborated in Sect. 2 can be seen in the Table 3. The accuracy increases as the number of input samples increases while the data outliers are filtered to improve the estimation. The number of input samples can be selected to reach a good tradeoff between the accuracy and computational speed. The actual values of the diagonal elements of stiffness matrix as the main stiffness parameters are listed in the Table 3. In the case that the exact values of eigenfrequencies and inertia parameters are accessed, the solution of Eq. (16) gives the exact values of stiffness parameters. The estimated values of stiffness parameters of the main diagonal of the matrix, by considering both the natural frequencies and inertia parameters estimation errors, are listed in the Table 3. Estimated M1 is the designation used to indicate that the first proposed stiffness estimation method which denotes the case that all the eigenfrequencies are accessed. Estimated M2 indicates that the second method of stiffness estimation is used, which is associated to the case that some of the frequency modes are not observable. In order to calculate Estimated M2 , two different cases in which only the first ten and eleven estimated modes are available are simulated, and the optimization problem defined by the LS estimator in Eq. (18) is solved. The results are listed in the Table 4. The estimation error of the associated underdetermined LS estimator reduces as more frequency modes are known. It is worth noting that our extensive simula- tions show that high values of error in the input data of stiffness estimation problem, namely the estimated eigenfrequencies and inertia parameters, may cause instability of the nonlinear numerical solver in Eq. (16). From the stability perspective, the LS estimator outperforms the first method which is based on solving the nonlinear Eq. (16), because the eigenfrequencies of Λ in the LS estimator are forced to be positive, which ensures the convergence of the solver. In other words, Estimated M2 is more robust to the input data, and is recommended also when all the frequency modes are available.

Loads, contact stress and accumulated damage
Contact load and stress validation on the three different gear stages The input and output torques derived in Sect. 2.5 and shown in Eq. (25) can be seen in Fig. 6. There one can see that all torques present similar oscillating pattern around their nominal values, and that there is little difference between the output and input torques at adjacent stages, see Fig. 6d-b and e-c. One can estimate the stresses at a sun-planet gear pair by inserting the input torques shown in Fig. 4 at the first order approximation using Taylor series, based on Eq. (12). The results for such operation can be seen in Fig. 7, together with results from high-fidelity Simpack multi-body simulation model of NREL 5 MW drivetrain reported in [12]. The results show reasonable agreement between the results for the simplified torsional model and the high-fidelity multibody simulation platform. (Fig. 7).

Calculation of average accumulated damage
The SN curve parameters for pitting fatigue damage calculations depend on material, load and gear geometry. These two parameters for the pitting fatigue damage estimated in this paper are obtained based on ISO 6336-2:2019 [14]. Only one region is assumed for the SN curve. For instance, k is 3.051 × 10 44 , and m is 12 for the sun gear of the 2nd stage. The numbers of stress cycles of sun gear is multiplied by 3, because it meshes with 3 gears simultaneously at each revolution. The accumulated damage of 1st stage sun, 2nd stage sun and 3rd stage pinion gears due to contact fatigue for one hour of operation by using the proposed digital twin-based remaining useful lifetime monitoring approach is listed in the Table 5. The online accumulated damage D of the 2nd stage sun gear is shown in the Fig. 8.

Practical implementation challenges and realtime performance
Challenges to overcome for practical implementation In practical turbines, uncertainties are an important issue. The different sources of uncertainties and their influences on the digital twin-based degradation estimation approach grounded on drivetrain torsional equivalent models for estimation of residual life of the drivetrain shafts are investigated by Moghadam and Nejad in [24]. The latter is performed by using stochastic models and statistical approaches, where uncertainties in both the model estimation and real-time measurements are addressed to mitigate their influence and improve the efficacy of the digital twin approach in presence of uncertainties. The statistical uncertainties in fatigue calculation due to material uncertainties can be accounted for by stochastic modelling of damage as explained by [24]. It is also worth noting that the employed least-square estimator is robust to the measurement noises, so that it can mitigate the influence of measurement noise in general, and cancel the influence of Gaussian noise [25]. The other issue in practical turbines is a call for computationally fast and inexpensive health monitoring techniques. The proposed algorithm is based on inexpensive linear equivalent models of the drivetrain, and the computational time depends on the sampling frequency and the number of input measurement samples. The algorithm can be managed to be executed only in order of a second. The method relies on a linear regression-based estimator and a few measurement samples, linear torque observers, and the real-time cycle counting of the estimated stress. The selected linear 14-DOF model of drivetrain is a compromise between simplicity and accuracy, which demonstrates similar results to results of high fidelity models for estimation of load aimed at estimating the degradation of the gears. The practical implementation of the algorithm integrated with turbine fully automated control and monitoring systems means the automated pre-processing of the turbine measurements to continuously feed updated input data into the proposed system identification approach to estimate real-time equivalent model of drivetrain in the digital twin framework. This task aims at the optimization of data streaming and continuous processing architectures to deal with the experimental data aspect of the digital twin framework. Within the preprocessing stage, a significant step will be on time signal correlation methods to verify and ensure time synchronicity when the data comes from different sources by mapping them on common spaces and dealing with the effect of different samples rates. Another practical issue with the realization of this idea is the possible need for new sensor installations to access additional torsional measurements, communication links with the main operation control unit, and the processing power for the execution of algorithm in near real-time.

Real-time performance metrics
As discussed in detail in Sect. 2, the digital twin algorithm consists of three main components, namely drivetrain equivalent model identification, load estimation and damage calculation. Therefore, for the simulation-based studies of this paper, the run-time of algorithm can be estimated by the summation of data recording time, time for estimation of equivalent model parameters, time for estimation of load and time for calculation of damage. The model identification computation time, t model ; is mainly influenced by the linear regression estimator operation, which its computational time depends on the degree of model complexity and the number of input samples. For the case of 14-DOF model, the model identification algorithm can estimate the model parameters by less than 5% error only by using 14 data samples, where the estimation error can be reduced by increasing the number of input samples. By assuming 30 data samples, with the sampling frequency 300 Hz, the time length of recorded data block will be 0.1 sec. The processing time depends on the processing power, but can be managed to be only a fraction of a second. The estimation of load and stress computation time, t load ; is determined by the operation of designed torque observers which are simple arithmetic functions of estimated model parameters and measured real-time response, with the computational time very close to zero. The damage calculation computation time, t damage , is reliant on real-time cycle counting operation on the data blocks of stress time series, which can be executed quite fast in fraction of a second. Two real-time metrics are employed to evaluate the algorithm run-time, namely real-time factor (RTF) and run time (TR), as they can be estimated by [26] TF = t process t input ; TR = t input + t process ; t process = t model + t load + t damage : These equations can be used to provide a rough estimation of the real-time performance of the algorithm, based on testing it with a system with 1.9 GHz Intel Core i7 CPU and 16 GB memory, which leads to RTF = 0.8 s + 0.1s + 0.2s 0.1 s = 11; TR = 0.1s + 0.8s + 0.1s + 0.2s = 1.2s: Generally speaking, for a system to be considered realtime, RTF should be ≤ 1 [27]. The run-time of the algorithm can be compensated in the proposed digital twin approach by repeating the results of data processing phase for all the data blocks captured during the processing time. The latter is expected to provide fair enough results for the drivetrain components' fatigue damage estimation, since degradation usually happens during longer periods of time than one second. Hardware-in-the-loop (HIL) simulation of the algorithm to study the possibility of executing the algorithm in real-time in practice is considered as the future work.

Conclusion
The possibility of using 14-DOF lumped parameter model as the drivetrain digital twin model for monitoring the remaining useful lifetime of the gears in the different gearbox gear stages due to contact fatigue stress was investigated. An algorithm for near real-time estimation of the parameters of drivetrain equivalent model by using the torsional measurements was proposed and tested by simulation studies. The estimated contact load and stress obtained by using the designed load observers in the proposed digital twin approach were validated by the results obtained from a Simpack high-fidelity simulation model, which showed a fair agreement between the results, whereas the proposed digital twin approach based on a linear torsional model is computationally fast and can be implemented integrated with turbine fully automated control and monitoring systems for online monitoring of gearbox components. The estimated stress was later used for near real-time estimation of the fatigue damage of the gears by using a physics-based degradation model. The influence of pitting faults on the system's resonances was demonstrated by simulation studies, which can support the fault diagnosis of the gearbox.
Hardware-in-the-loop (HIL) simulation of the algorithm to check the possibility of executing it in real-time for the fault prognosis purpose, tackling with the described practicality issues and implementation of the algorithm in an operational wind turbine drivetrain system are considered as the next steps that will be investigated in the future work.
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.