Development of a digital twin of an onshore wind turbine using monitoring data

The present paper has two main objectives: (1) to create a digital twin of an onshore wind turbine to provide to the owner of the wind farm a tool for continuous tracking of accumulated fatigue damage and evaluation of alternative operation strategies; (2) to perform the first tasks for the creation a reliable numerical of a floating wind turbine to simulate experimental data for the testing of monitoring tools on this type of wind turbines. A numerical model of an onshore wind turbine was created with FAST software developed by the US National Renewable Energy Laboratory (NREL). The structural and mechanical properties of the structure were calculated from the geometric properties of the tower and blades. The aerodynamic properties of different sections of the blades were computed with 2D models created with ANSYS Fluent. Having the wind turbine fully characterised, the control mechanisms were calibrated using data from the manufacturer catalogues. Numerous time series of wind excitation for different operation conditions were generated, and the structural response computed. Additionally, the Supervisory Control and Data Acquisition (SCADA) operating in parallel with a series of extensometers and accelerometers located in strategic points along a real wind turbine tower and blades provided a full description of the internal forces in both locations and the estimation of the modal properties for different operating conditions. Both measured and simulated response allowed the identification and validation of structural dynamic properties and static and dynamic internal loads. In spite of several approximations and simplifications, the results obtained with the numerical model were in good agreement with the ones measured in the field.


Introduction
The main objective of this work is to set the way to a reliable off-shore wind turbine numerical model which can then be used to simulate the structure response and test alternative monitoring strategies algorithms for data processing, overcoming the lack of experimental data available for this type of structures. In a first approach, it is fundamental to develop a model that can adequately reproduce the behaviour of onshore wind turbines, where results can be validated with experimental data. For that purpose, a wind turbine located in Tocha, Portugal, monitored in the context of the WindFarmSHM project with the sensors described in detail in [1] was chosen as case study. This paper describes how the numerical model was obtained as well as the different validations made by comparing numerical results with experimental ones.

Numerical Model
The modelled wind turbine is a 100 [m] rotor diameter and 93. 3 [m] tower height VESTAS V100 with 1.8MW rated power capable of working between 4 [m/s] and 20 [m/s]. The accuracy of the results obtained is strongly dependent of both structural properties and loads acting on the structure. While the first point can be obtained without great difficulty the former is mainly influenced by the aerodynamic properties of the blades. To accurately model both effects it was used FAST software, therefore a brief description of the software main principles is included.

Structural Properties
The static and dynamic response of the structure to any applied load are mainly dependent on the behaviour of the tower and blades. As FAST uses a generalized degree of freedom approach based on the modal configuration of the tower, a finite element model was made to obtain these configurations.
In this initial FEM model, the tower was simulated providing the linear mass density and SS/FA stiffness in 37 points distributed along its height. The blades were discretized in 49 different sections and to each section were assigned 2 directions stiffnesses, density and structural twist. As only the external geometry of the blades and its global mass were known, the stiffness properties and variation of linear weight were scaled from the 5MW baseline turbine available in FAST certification tests [2]. The chord scale (k C ) was calculated as the ratio between the chord of the VESTAS V100 and the reference blade. Assuming these can be treated as thin wall hollow sections then the mass scale can be obtained directly from the product of the chord and thickness scale (k t ). Matching the global mass obtained for both, allowed to calculate the latter. With this assumption, the inertia scale (k I ) for each section can finally be calculated as the product of the 3 rd power of the chord scale by the global thickness scale. The foundation was modelled with shell elements with properties taken from final design drawings. Based on these information, a first FE model was built in ROBOT Structural Analysis (RSA) assuming a rigid connection between the nacelle and the hub. As the eigenfrequencies of this model did not match the experimental ones, two main sources of error were identified, being the blades stiffness and the drivetrain flexibility. While the first can be adjusted by changing the blades material properties (equivalent to an additional global inertia scale factor), the latter can be introduced in the model by means of an elastic connection between the rotor and the hub.
Having the dynamic behaviour of the structure tuned, the two principal modes configuration for the tower in each direction were deducted from the tangent line at the tower base so that they respect FAST assumption of 0 slope at this section and fitted to a 6 th degree polynomial with unitary displacement at the top of the tower, as required [3]. In FAST model the foundation was assumed to have infinite stiffness to translations and rotation about the tower yaw axis and the remaining component of the stiffness matrix were evaluated from the bending moment and rotation of tower base of the described FEM model. To simulate this effect, it was used the HydroDyn module with an imposed linear stiffness matrix only.

Aerodynamic Properties
The aerodynamic properties of the blades are an input of FAST and can be calculated with carefully built CFD models. The CFD model properties used to obtain the lift and drag forces (F L and F D respectively) and corresponding coefficients for the different angles of attack (α) for the adopted aerofoils were previously validated with the well-documented NACA0012 aerofoil, which was study with wind tunnel experiments [4]. In order to capture the variation of the aerodynamic properties along the blade, it was divided in 10 different sections equally spaced - Figure 1 -and a CFD model was made for each, for different angles of attack. The adopted conventions are defined in Figure 2.  Figure 1.
Representation of the adopted blade sections. Figure 2. Representation of the adopted forces conventions for a given wind (U 0 ) and rotor speed (Ω) [5].

CFD model validation
The CFD simulation was carried out with ANSYS/fluent code based on the Reynolds Average Navier Stokes Equations (RANS). As described in [6] the Spalart-Allmaras (SA) turbulence model [7] has been systematically tested over the years and provided accurate results, at least for the precision assumed necessary in this work. It was also used the SST k-ω turbulence model, as it captures both the accurate formulation of standard k-ω model in the near-wall region and freestream independence of the k-model in the far field [8] . The aerofoil was treated as a no-slip wall and all the equations were solve using second order schemes only. A preliminary convergence analysis was made using 5 different model orders, where the height of the computational domain presented in Figure 3 was divided in 200 · i elements, being i the model order, having obtained the drag forces presented in Figure  4. Although mesh independence appears to be obtained for the 3 rd order model, at least using SST, as computations for higher angles of attack may require a more refined mesh, the 4 th order model was chosen. This is represented in Figure 5.   Reynolds and Mach numbers similarity with the experimental results obtained by Ladson in [9]. Figure 6 represents these coefficients for all the tested angles of attack. Both models are in good agreement with the experimental results used for comparison and while the lift coefficient before stall is mainly Reynolds number independent the same does not apply to the drag coefficient where higher Reynolds number lead to lower drag forces.

CFD results
As the simulations for angles of attack (AoA) higher that stall angle become more complex it is important to estimate the values that the different sections will experiment during the wind turbine operation conditions. This evaluation can be made using Blade Element Theory (BEM) expressions -equations 1 and 2 (see notation in Figure 2).
However, as the forces decomposition depend on the axial and tangential induction coefficients (a and a respectively) that by turn depend on the forces, an iterative process is needed. Defining the blade solidity (σ r ) as the ratio between the blades area and the rotor area it is possible to write an expression for both induction coefficients based on the parameters of the previous iteration -equations 3 and 4, where C and C ⊥ represent the component of the lift and drag forces in and out of the rotor plane respectively.
In order to evaluate the angles of attack that will be necessary for each blade, it was assumed a linear behaviour near the 0º angle of attack so that the lift coefficient could be computed based on equation 5 with slope (m) and zero lift angle of attack (α 0 ) obtained from CFD simulations between -5º and 5º. It was also assumed the expression for the drag coefficient obtained for a fully laminar boundary layer over a flat plate so that its value depends on the Reynolds number only and can be calculated with equation 6 [10].
Knowing the structural twist distribution and choosing the pitch angle in accordance to the SCADA system records, it is possible to obtain the angles of attack for the 9 last sections presented in Figure 7 1 . The expected angles of attack are within a rather short range which can be analysed with the CFD model. All sections were tested for AoA between -10º and 12º and Reynolds number equal to 2 × 10 6 , which was chosen as a representative value for the different operating conditions. Having the drag and lift coefficients for these angles of attack two important corrections should be considered. The first one is to correct the values to take into account the 3D stall delay caused by a rotating blade. This can be achieved by the methods proposed by Selig in 1998 [11] and Eggers in 2003 [12] for the lift and drag coefficients, respectively. Furthermore, as only angles of attack in a restricted domain were tested and a full description for any other value 1 The first section is not analysed as it was considered a circular section with constant lift and drag coefficients.  [13]. The representation of lift and drag coefficients for one of the aerofoils can be seen in Figure 8.

Control Systems Simulation
The control system properties were defined using the simplified routines presented in [14]. As there was no available information about the pitch angles below rated speed, it was necessary to look into the SCADA data and define an average minimum pitch angle that was fed into the pitch control routine beside the rated rotor speed. The torque control was defined so that the power output would match the theoretical values, setting the appropriate values of gear box efficiency, power law coefficient of Region 2 and Rated Torque value of Region 3 2 . As Region 2 is assumed to have constant pitch angle and constant speed, an infinite number of equilibrium values between aerodynamic torque and resistant torque can be defined associated with different rotor speed and hence different thrust forces due to different angles of attack. Having this in mind, the power law coefficient was defined by an iterative process in order to minimize the differences between theoretical and numerical thrust forces. The gear box efficiency can then be 2 The mentioned regions are in accordance with the nomenclature used in [14]. Region's 2 torque is assumed to be proportional to the 2 nd power of generator speed and Region's 3 torque is assumed to be constant. calibrated matching the theoretical power output with the generator speed obtained in the last iteration of the previous step. The rated torque value can be defined from the already determined gear box efficiency and rated speed so that the theoretical rated power is achieved for wind speed above 12 [m/s]. The power and thrust curves obtained from 10 minutes simulations under nonturbulent steady wind conditions 3 for different wind speeds can be seen in Figure 9, and a good match between numerical and theoretical data is obtained.

Results
With the model fully characterised, it is now possible to compare the numerical and experimental data available. In order to validate the simulated response, a comparison of some characteristic curves can be made with the SCADA data and further validations can be obtained comparing the internal forces and dynamic properties derived from simulated and measured datasets.

Single realizations comparison
To check that the modelled dynamic behaviour is the pretended, both experimental and numerical natural frequencies can be compared using operation modal analysis techniques. The response of the structure was simulated during 30 minutes with the rotor parked and then divided in shorter series of 300 [s]. The power spectrum was calculated for each 5 minutes dataset separately and then averaged. The spectrum obtained from both experimental and simulated accelerations at 75% of the tower height can be seen in Figure 10, where it is marked a reference value for the resonant frequencies obtained from experimental data. The non-calibrated model is also represented for comparison. While the 1 st mode could be accurately predicted from the original model, the 2 nd mode presents bigger deviations from experimental data due to its strong dependence on the interaction between tower and blades, only accounted for in the calibrated model. However, and even in this model, although the principal resonant frequencies are correctly identified, other peaks of smaller amplitude are not correctly obtained. Based on the experimental data and the FE model, these frequencies are associated with mode shapes mainly dependent on the blades behaviour and so these less accurate results can be explained by all the simplifications made in the structural and mechanical properties of these elements.  The numerical power, rotor speed and pitch curves can be superimposed with the ones recorded by the SCADA system as represented in Figure 11. Analysing the rotor speed and pitch angle curve it is clear that the control implemented reproduced adequately the values 8 obtained from SCADA system for most wind speeds. The results associated with lower wind speeds are not a major concern, as they are associated with low forces 4 . The power output seems to be in good agreement with the expected value for all the tested scenarios. Finally, and to compare the internal forces of the model with the ones measured in the field, a selection of some 10 minutes interval with wind speed oscillating around an approximately constant mean value was used. To select such intervals, the SCADA data from February 2019 was analysed and 10-minute time intervals with the properties presented in Table 1 were selected. Based on these wind characteristics it was possible to construct 3-dimensional wind grids with equivalent mean speed and turbulence intensity with TurbSim pre-processor from NREL [15].  Table 1.

Extended results
To provide a more robust validation of the model, 100 independent 10 minutes simulations, considering a 10% turbulent intensity, for different mean wind speed were computed. The thrust force acting on the rotor, and corresponding bending moment at the control section, obtained for each computation with a given prescribed mean wind speed is plotted in Figure 13, where the error bars identify the maximum and minimum mean values obtained, considering a 0.5 and 1 tower drag coefficient (a limiting value taken from [16]). As can easily be checked, the lower drag coefficient leads to a better agreement with experimental data.  Figure 13. Comparison between numerical and experimental forces for 100 independent generated time series.
In the same way, the response power spectrum was computed for each independent data set and then averaged for a given wind speed, resulting in the 3D array represented in the colour maps of Figure 14. From the analysis of the stability diagram for a fixed mean wind speed it can be seen that the principal frequencies are stable for all simulations and while some of the higher values can easily be matched with the ones from Figure 10, it is clear that some previously absent frequencies are now important for the response content.   10 The nature of this peaks can be easier understood when the frequencies evolution with the mean wind speed is analysed. In these diagrams, the natural frequencies of the tower can be seen to be stable over all the simulations but the newly identified peaks do increase with the mean wind speed until a given value for values near the rated wind speed. In fact, this peak is due to each blade passing through the non-uniform turbulence field and hence its frequency of 3 times the rotor speed (0.745 [Hz] for the 14.9 [rpm] associated with rated wind speed).

Conclusions
The developed model was able to make suitable predictions of the thrust force and power output using only data from project drawings and theoretical curves. To obtain good results regarding the dynamic behaviour of the structure it was fundamental the calibration with the experimental data, which allowed to correctly identify the principal global modes of the structure. The frequency peaks associated with blades local modes were not identified correctly but the global response is in good agreement even with the simplifications and assumptions that were needed. In future works it is important to define the blades with more correct information and extend the validations made to local internal forces and modal response of the blades. Finally, the bending moments at the tower predicted by the model were higher than the ones measured on the field for the transition region between torque and pitch control, which can result from either a too simplified control algorithm or an overestimated wind tower loads, so future studies are needed on this matter. Considering the actual interest of wind farm owners on the continuous evaluation of accumulated fatigue damage and evaluation of alternative operation strategies through the use of a digital twin fed with SCADA data, this modelling exercise is also being extended to fulfil this need.