Aeroelastic validation and Bayesian updating of a downwind wind turbine

Downwind wind turbine blades are subjected to tower wake forcing at every rotation, which can lead to structural fatigue. Accurate characterisation of the unsteady aeroelastic forces in the blade design phase requires detailed representation of the aerodynamics, leading to computationally expensive simulation codes, which lead to intractable uncertainty analysis and Bayesian updating. In this paper, a framework is developed to tackle this problem. Full, detailed aeroelastic model of an experimental wind turbine system based on 3-D Reynolds-averaged Navier-Stokes is developed, considering all structural components including nacelle and tower. This model is validated against experimental measurements of rotating blades, and a detailed aeroelastic characterisation is presented. Aerodynamic forces from prescribed forced-motion simulations are used to train a time-domain autoregressive with exogenous input (ARX) model with a localised forcing term, which provides accurate and cheap aeroelastic forces. Employing ARX, prior uncertainties in the structural and rotational parameters of the wind turbine are introduced and propagated to obtain probabilistic estimates of the aeroelastic characteristics. Finally, the experimental validation data are used in a Bayesian framework to update the structural and rotational parameters of the system and thereby reduce uncertainty in the aeroelastic characteristics.

and significant differences with respect to experimental results. In this regard, we explore the use of robust reduced-order models (ROMs) for wind turbine applications to balance the computational overhead and maintaining acceptable accuracy.
In the academic community, three-dimensional (3-D) Navier-Stokes aerodynamic/aeroelastic simulations for wind turbines with structural details of the wind turbine such as tower and nacelle have been considered in only a few studies. [5][6][7][8][9] Multiple investigations have been performed considering only the aerodynamics, [10][11][12][13] while structural analysis of rotor blades with complex geometry and material compositions has also been considered. [14][15][16] The aerodynamic computations have generally been performed with some limitations, for example, with a simplified structural representation. One of the first 3-D simulations of wind turbine rotors at full scale was reported with an isogeometric analysis based on non-uniform rational B-splines, 17 employed for the fluid-structure interaction (FSI) modelling of the NREL 5-MW offshore rotor. However, the effect of the tower and nacelle was ignored. This study was followed up with a more detailed representation of the rotor-tower interaction 18 and then demonstrated in a Windspire vertical axis wind turbine (VAWT). 8 The type of coupling between the computational fluid dynamics (CFD) and structural solvers in aeroelastic simulations determines the consistency in forces and accelerations. In this regard, there have been developments based on loose coupling 19 and a tightly coupled method. 4 However, these simulations consider upwind configuration of the wind turbine. In this paper, we consider a downwind configuration that has to deal with more complex flow features since the blades are subjected to the tower wake at each rotation.
Another issue affecting the accuracy of wind turbine simulations is related to the inherent uncertainties in the structural parameters and the imposed fluid boundary conditions that have to be considered in the numerical model. One of the first probabilistic treatment of parameters in wind energy was performed by estimating the uncertainty in wind turbine power output and annual energy production using the Weibull distribution. 20,21 The uncertainty in wind direction was shown to be a main source of discrepancy in assessment of wind farm power output obtained from numerical simulations, when compared with experimental data sets. 22 Monte Carlo-based approaches have also been used to sample the effect of multiple parameters on the wind power output. [23][24][25] A reliability-based design optimisation of wind turbine blades under wind load uncertainty was carried out based on 249 groups of wind data 26 to consider the variation in wind. Common UQ propagation techniques, such as polynomial chaos expansions and stochastic collocation methods, have also been used to propagate uncertainties in wind turbine applications. 27,28 Most of the probabilistic studies in the field of wind turbine simulations have been based on simplified aerodynamics. More recently, 29 an aerodynamic shape optimisation of wind turbine blades using a Reynolds-averaged Navier-Stokes (RANS)-based fluid model and an adjoint method was presented. However, unsteady effects were not considered due to computational limitations. An aerodynamic design optimisation for a horizontal axis wind turbine (HAWT) under geometric uncertainty was studied 30 using the univariate reduced quadrature (UEQ) approach, but the aerodynamics was based on blade element momentum (BEM) and a single airfoil shape was employed for the entire length of the blade.
A stochastic analysis of flow-induced instabilities because of uncertainties in flow forces as well as structural properties was performed 31 using a linear stability analysis and modelling the aerodynamics based on Theodorsen theory. The effect of this randomness on the onset of instability is clearly observed, including observance of coupled-mode flutter at speeds below the designed operational speed. Considering the importance of including uncertainties in numerical simulations, in this paper, the input parameters in the wind turbine are treated as uncertain.
As already mentioned, the computational expense of full-order aeroelastic models makes uncertainty propagation intractable. As such, data-driven techniques are being increasingly explored in order to obtain reliable estimates at low costs. Calibration of simpler models (RANS, Jensen wake models) using high-fidelity large eddy simulation (LES) data has been implemented to reduce the computational expense. 32,33 The model parameters of a simplified finite element-based structural model are calibrated using Bayesian inference in order to predict the blade dynamics. 34 Modal decomposition methods have also been applied to high-fidelity simulation data in order to obtain simplified models for studying wake dynamics. 35 Dynamic mode decomposition has been used to build a reduced-order model (ROM), 36 which was then embedded in a Kalman filter to produce a time-marching algorithm. This dynamic model could also use new data to dynamically update the ROM to provide real-time estimates. Building on data-driven approaches, in this paper, we explore the development of aeroelastic ROMs based on training data from high-fidelity solvers and employ them to propagate uncertainties.
The contribution of this paper is threefold. First, dynamic numerical simulation of a full-scale rotating wind turbine of downwind orientation is presented, considering detailed structural components such as blade, tower, and nacelle with a high-fidelity fluid solver and detailed aeroelastic characteristics are investigated. Second, a data-driven reduced-order modelling framework is developed in order to predict aeroelastic loads in wind turbines blades. Third, data from a wind turbine experiment are used to identify the structural and operational parameters of the wind turbine in a Bayesian framework. The novelty of this work is in terms of the detailed aeroelastic simulations for a downwind wind turbine, the development of ROM framework for wind turbine problems, and the use of Bayesian framework with experimental data in wind turbine problems. Accordingly, in this paper, we present a methodology to predict dynamic loading in a downwind-oriented wind turbine, considering uncertainties in both flow and structural parameters.
The framework is adopted from an earlier study by the present authors, 37 applied to the case of aircraft wing. Initially, a high-fidelity aeroelastic solver is developed with a modal structural solver and RANS-based fluid solver. In terms of aeroelastic simulation of wind turbines, this is among the few simulations of a complete wind turbine of downwind orientation, considering all structural components. The aeroelastic solver is validated by findings of the wind turbine experiment, and accurate estimates of the blade moment variation with azimuthal position is provided by the solver. Extensive investigation of the flow-field in downwind wind turbines is presented. The aerodynamic solver, with prescribed forced motion of the blades, is used to generate data to train a system identification ROM for the fluid component only. The prescribed motion is generated through informative chirp signals that are expected to excite the frequencies of the blades and are important to predict the investigated aeroelastic phenomenon. The ROM predicts the dynamic aerodynamic loads on the wind turbine blades, which is then coupled with a structural solver. The coupled aeroelastic solver can be used to estimate the fluctuating loads leading to structural fatigue of the system. The dynamic aeroelastic characteristics are then investigated under blade structural uncertainty, which is a rapid process as the aerodynamic ROM does not need to be retrained. Finally, we use experimental data in a Bayesian inference setting in order to identify dominant source of uncertainty in a wind tunnel experiment.
The paper is structured as follows: Section 2 provides the theoretical description of the aeroelastic solver with details of the structural and fluid solvers, coupling scheme and interface treatment between stationary and moving frames of reference. Parametric investigations are performed in Section 3 with details about the effect of downwind configuration on the blade moments and loading. Section 4 is an overview of the ROM framework and details of the reconstruction of dynamic blade loads obtained with the ROM are provided. Uncertainty propagation is shown in Section 5 with use of data inference in order to identify stochastic parameters. Finally, Section 6 gives an overview of the findings of the paper along with recommendations for future research.

AEROELASTIC MODELLING FOR DOWNWIND WIND TURBINES
In this section, development and details of a high-fidelity aeroelastic solver are discussed. The wind turbine model is a three-bladed experimental wind turbine of downwind orientation developed by Verelst. 38 The experiments were performed at the Open Jet Facility of the Delft University of Technology. The rotor radius was 0.8 m with NREL S823 and S822 airfoils at inboard and outboard sections, respectively. There were two blade types: flexible and stiff. This is based on the number of layers of glass fibre used to construct the blades. Two and eight layers were used in the flexible and stiff set, respectively. The reader is referred to Verelst 38 for further details. Because of the internal glass fibre and foam-based structure, in addition to the twist and chord variation along the span of the blade, a non-uniform distribution of mass and stiffness properties must be considered in the structural model.

Development of structural solver
The structural model of the blade must take into account the non-uniform structural properties s of the blade, which will be treated as uncertain in the later part of this paper. A modal structural solver is developed using an ordinary differential equation (ODE) model of the wind turbine blade undergoing rotation. Approximate sectional mass and stiffness properties are obtained from the experimental measurements, which were based on an optimisation procedure with the measured eigenfrequencies as objective function. 38 The equation of motion of a rotating beam, FIGURE 1 Schematic of the three-bladed downwind wind turbine with nacelle and tower of Verelst 38 with coordinate system shown in Figure 1, and with non-uniform sectional properties is given by 39 : where u and v are the edgewise and flapwise displacements in x and y respectively, while is the torsional deflection about z-axis. E and G are the Young modulus and shear modulus, I x and I y are the moments of inertia, J is the polar moment of area, r x and r y are the distances between centre of flexure and centre of gravity in x and y, m is the mass per unit length of the beam, while R is the distance from centre of rotation to the root of the blade and Ω is the rotational speed of blade. The reader is referred to Carnegie William 40 for further details.
Flapwise, edgewise, and torsional modal frequencies and mode shapes are obtained from (1 to 3) for the wind turbine blades. The frequencies of vibration are validated with the experimentally obtained natural frequencies; see Table 1. It can be seen that the first ( 1 ) and second frequencies ( 2 ) are close to the experimental values. It should be mentioned here that in the experiments, edgewise modes are not reported.
For the flexible blades, the third frequency 3 does not correspond to the experimental observations. This is attributed to the fact that the third frequency in the simulation is from an edgewise mode. However, note that 4 is a good match with the experimental value for the third frequency, which is shown as 4 here. Based on this, the structural model can be considered as validated and the estimated modes are used to build the modal structural solver.
We express the structural equations (1)(2)(3) in the modal form to reduce dimensionality of the model, by assuming the solution to be given by the following: where U n (z) is the n th mode shape obtained from modal decomposition and f n (t) is the generalised displacement. Corresponding to the frequencies obtained in Table 1, four mode shapes are considered in this investigation. Furthermore, since we are interested in the aeroelastics, we show the Full mesh of the computational domain: only rotating part results with a flexible blade (Flex 1) for the rest of the paper. Substituting (4) in Equations (1) to (3) and after simplifications (not shown here for conciseness), the modal structural equation reduces to the following: where Q n (t) is the aerodynamic forcing term, to be interpolated from the fluid to the structural mesh.

Mesh generation and fluid solver
We use a standard hybrid mesh for the present investigation, with a structured mesh for resolving the boundary layer and an unstructured mesh for the rest of the domain as shown in Figures 2 and 3. Such hybrid meshes are known to be able to provide accurate results for wind turbine flows. 41,42 The unsteady compressible Navier-Stokes equations are discretised using a cell-centered, finite-volume flow solver. A second-order upwind scheme 43 is used for spatial discretisation, and temporal terms are discretised with a second-order implicit time-stepping scheme with dual-time sub-iteration. The solver employs a low-Mach number time-accurate unsteady preconditioner, where physical time derivatives march the unsteady terms, while the preconditioned time derivatives are used for numerical discretisation and iterative solution. 44,45 We use a RANS-based turbulence model -SST k − , 46 which has been previously employed for wind turbine applications. 4,47,48 The boundary layer is resolved on the walls of the blades, nacelle, tower, and all the blade attachments, ensuring y + < 1. The details of the mesh and computational resources are mentioned in Table 2, where the time step Δt is stated in terms of the rotation angle of the wind turbine blades.
Despite the use of RANS-based closures in earlier wind turbine investigations, it should be mentioned here that higher fidelity methods such as LES and direct numerical simulation (DNS) are expected to provide more accurate representation of the loads. However, the computational costs of LES and DNS are prohibitively expensive to be employed in the ROM framework employed in this paper. The limitations of the RANS-based models to match experimental measurements have been reported in the MexNext report, 49 where CFD estimates from different solvers have been compared with predict the dynamics of the MEXICO rotor. 49 Additionally, a summary of CFD solver details for different investigations of the MEXICO blade is provided in Carrión et al. 50 In the context of the current study, a comparison of performance of turbulence models is not provided here. Since the SST k − model has been employed in earlier investigations, this model is chosen here. With the development in computational power, this presented framework can be extended to high-fidelity fluid models, such as LES.
For the rotating domain, the Navier-Stokes equations are written in the moving frame of reference in the relative velocity formulation, with addition of a Coriolis and a centrifugal force term. 51 Because of presence of the tower, the computational domain consists of a rotating and a nonrotating domain. A sliding mesh technique is employed for the interfacial treatment between the two domains. In the aeroelastic setting, the fluid and structural meshes are nonconforming and the structure is 1-D. We employ infinite plate splines (IPS) 52 to map the structural mesh  Computational domain in millions of cells and CPU hours consumed per rotation of the wind turbine displacement to the fluid mesh. Using the principle of virtual work for conservation of energy, the forces from the fluid mesh are mapped onto the structural mesh Q n (t) in (5). In order to avoid instability because of sudden forcing at the beginning of the simulation in (5), an initial damping term is introduced (for about quarter rotation of the turbine) through a damping coefficient, which is gradually reduced to zero.

VALIDATION OF WIND TURBINE DYNAMICS BY THE AEROELASTIC CODE
The experiment of the downwind system 38 was performed at varying operating conditions: wind speed, Tip Speed Ratio (TSR) of the turbine, temperature, and pressure of incoming wind. Additionally for each experiment, the revolutions per minute (RPM) was fluctuating slightly for each rotation. In order to obtain better match between experiments and aeroelastic approximations, we select data with minimal variation in rotation speed, which remains approximately 3% for the investigated cases. Motivated by this, later, this rotation parameter will be treated as uncertain.
In this section, we provide a detailed analysis of the performance of the aeroelastic solver in estimating the flow physics for a downwind turbine. Parametric effects on the wind turbine performance are also assessed and agreement of the aeroelastic code with experimental data is established. In terms of computational expense, approximately 362 CPU hours are consumed for 1 rotor rotation of the wind turbine, which is considerably high for performing UQ. In the next section, we develop ROMs to predict wind turbine aeroelastics.

Agreement with experiment
In the experiment, 38 strain gauges were fixed to the glass fibre-reinforced beam on the blade at the root and at 30% radial position. They were calibrated by statically loading the beam at different radial positions, and a linear function relating the output of the strain gauge to the bending moment in the blade was defined. In this regard, we can represent the blade moment in two forms: (a) instantaneous aerodynamic moment derived from the blade pressure distribution and (b) structural moment derived from blade deformation, which is smoother due to the inertia of blade.
It can be said that (b) is filtered version of (a), where the forcing from the tower wake is represented with a smoother curve without large spikes. In the experiment, 38 the blade moments are reported in the form (b). Furthermore, significant reservations are reported with respect to the reliability of the absolute values of strain gauge measurements, whereas the moment trend with change in azimuthal angle is expected to be meaningful. 38 With this background, we present the validation of the aeroelastic code in terms of (b) above.
In Figure 4, the structural moment variation with azimuthal position is shown for two different experimental data sets at varying inlet wind velocity V and speed of rotation of the turbine N. Both the experiments were conducted at a low value of TSR, TSR ≃ 1.5 given by, where and R are the angular velocity and rotor radius, respectively. However, both were conducted at a fixed yaw angle and the flexible set of blades were used. The motivation for this choice is the simplification of the flow condition with a fixed yaw angle, and also, flexible blade measurements are more interesting for this research, which investigates aeroelastic effects. In both cases, the experimental structural moment measurements for each rotation of the rotor are plotted, ignoring the data from the first few rotations, and illustrated by the grey area in Figure 4.
It is clear that the spread in the measurements is considerable, with more than 50% variation for certain azimuthal positions. In particular, there is large spread in the moment when the blade is in front of the tower, which is around the 1.5 azimuthal position (marked by red line in figure), which is likely due to dependence on phase of shedding. However, a trend can clearly be observed for both the cases. The mean aeroelastic estimate follows the trend quite well, with domain frequencies and amplitudes being reproduced. In Figure 4A, the influence of the tower can clearly be observed; however, in Figure 4B, the effect is somewhat observed at a higher azimuthal position. This could well be related to the choice FIGURE 4 Comparison of structural moment estimated at 30% span location by aeroelastic code to experimental measurement. Red vertical lines indicate location of the tower in terms of the azimuthal position of the blade. TSR, Tip Speed Ratio of fluid model, which is based on RANS, that may not be adequate for predicting loads accurately. However, the large spread in experimental measurements shows that the uncertainty at these positions is high, and this ensures that there is some overlap with the computational estimates.
The source of this uncertainty will be clearer in the Bayesian identification problem, which will be discussed in Section 5.3.

Rotor and tower wake visualisation
The interaction of tower wake with blades leads to peaks in the generated moment (as seen in Figure 4), which could lead to structural fatigue.
This interaction leads to complex flow features as shown in Figures 5 to 7, where inlet flow velocity and speed of rotation are adjusted to obtain a TSR ≃ 6.3 in Figure 7A. The experiments with low rotational fluctuations were conducted at lower TSR values, which results in wind FIGURE 5 Vortex from blade tip, root, and nacelle of wind turbine at TSR ≃0.7

FIGURE 6
Vortex from blade tip, root, and nacelle of wind turbine at TSR ≃1.6. TSR, Tip Speed Ratio FIGURE 7 Vortex from blade tip, root and nacelle of wind turbine at TSR ≃6.3. TSR, Tip Speed Ratio turbine operating in stalled regime. At a higher TSR, the tip vortices are more compactly arranged, which can be captured more effectively due to higher mesh resolution in the near wake region. Figures 5A and 6A are vortex visualisations at TSR ≃ 0.7 and 1.6 respectively, which are the experimental test conditions. The isosurfaces are plotted with a Q-criterion 53 of 55.
For the higher TSR ≃ 6.3, the tip vortices are distinct and are preserved up to 1.5 rotations, after which the resolution in the mesh is not sufficient to capture the wake clearly further downstream. Aeroelastic influences are captured due to the higher resolutions in the near wake. For TSR ≃ 0.7, the tip vortices disintegrate quickly and cannot be clearly identified. This was also observed in Carrión et al and Li et al 4,54 and was attributed to the large pitch, which results in the vortex reaching the coarser mesh region earlier. However, in the present investigation, the wind turbine is of downwind orientation, in which case the interaction of the tower wake with tip vortices is even more significant. For TSR ≃ 6.3, it can be seen that the tower wake influences the disintegration of tip vortices in the lower half of the blade sweep area (below the nacelle), resulting in weakening of the vortices in that region. The tower wake growth is not entirely evident in this region at this value of Q and is attributed to the lower inlet velocity. As a result, tower vortices are much weaker than tip vortices.
By contrast, the tower wake in case of Figure 5A can be clearly observed. In Figure 6A, the inlet flow velocity is close to that in case of TSR ≃ 0.7. It can be seen that the tower wake in Figure 6A is prominent and is in agreement with Figure 5A. It can be concluded from these observations that the tower vortices are much more influential at lower TSRs. Another interesting observation from the vortex visualisation in Figure 7A is the production of large-scale vortices from multiple locations in the mid-span of the blade. This was also observed in Zhang et al 55 and is attributed to the large change in angle of attack in the trailing edge at that location. Also, there are large vortices developed from the nacelle, and a strong interaction with the root vortex in the downstream of the axis of the turbine is observed.
The difference in the growth of the vortices in the three test conditions can clearly be observed in Figures 5B to 7B, where the vorticity magnitude contour is plotted in a plane parallel to the flow direction and bisecting the tower along its length. For TSR ≃ 6.3, the shedding of the tip vortices can be observed distinctly for approximately 2 rotor rotations, after which it disintegrates. The root vortex, nacelle, and mid-span blade shedding are also clearly visible. Close to the tower, a strong circulation is observed, and the interaction of tower wake with the tip and root vortices can be clearly seen. In Figure 5B, the lowest TSR ≃ 0.7 is plotted, where stronger tower vortices are observed, and as observed with the 3-D visualisation, the tip vortices are not distinct. Similarly, for TSR ≃ 1.6, strong tower vortices are clearly observed, with a comparable inlet flow velocity. However, in this case, it is interesting to note that the tip vortex shedding can be seen for about 1 rotation. Because of the higher rotational speed of the blade, the pitch of the tip vortices is smaller compared with TSR ≃ 0.7, and the mesh resolution is still fine enough to capture the wake.

Low frequency variations in blade moment
For further investigations, we consider TSR ≃ 1.6 and 6.3 to compare the flow field characteristics, since they demonstrated significantly different wake regions. In Figure 8, the aerodynamic moment at 30% span location is plotted for 12 rotations of the turbine at these two TSR values.
The peaks in the moment history correspond to the instants when the blade is in front of the tower. At TSR ≃ 6.3, the blade is subjected to slightly higher mean moment (that is the moment when it is not in front of the tower). An interesting observation is that there is a low-frequency unsteadiness in the moment peaks for both cases. For the lower TSR ≃ 1.6, the amplitude of this low-frequency unsteadiness is observed to be significantly higher with a lower frequency of about 0.37 Hz compared with 1.22 Hz at TSR ≃ 6.3.

FIGURE 8
Comparison of aerodynamic moment histories at TSR ≃ 1.6 and 6.3. TSR, Tip Speed Ratio Spontaneous low-frequency modes are not uncommon in fluid dynamics, even in simple geometries, for example, in flat plate 56 and circular cylinder, 57 where they are attributed to the shrinkage and enlargement of the recirculation region. Both of these studies classify the wake formation temporal history into a high-energy mode H and a low-energy mode L, which are characterised by higher and weaker fluctuations in the shear layer respectively. In Najjar and Balachandar, 56 the modulation in the amplitude is associated with the compact spanwise vortices close to the wall, resulting in shorter recirculation region. This can be associated to the difference observed in the low-frequency characteristics at the two TSRs. Stronger vortex strength is observed at TSR ≃ 1.6 compared with TSR ≃ 6.3, which may be attributed to the formation of an analogous H regime as in Najjar and Balachandar. 56 This complex phenomena can also be influenced by the rotation speed of the turbine, as the low frequency in both the TSR cases are also different. Since the objective of this investigation is not to characterise this low-frequency behaviour, it is not pursued further here.

Effect of blade elasticity
The change in the blade forcing because of inclusion of blade aeroelasticity is interesting since the nature of the forcing could lead to structural fatigue. Multiple simulations are carried out considering rigid, rotating blades (referred as rigid) and then with flexible rotating blade (referred as aeroelastic). The simulations are allowed to run for around 5 rotor rotations in order to eliminate effects of the initial condition. In Figure 9, comparison of aerodynamic moments obtained at 30% span location with rigid and aeroelastic codes is shown for TSR ≃ 1.6 and 6.3. The aeroelastic effect slightly varies for each rotation cycle, and hence, a dominant pattern is not explicitly observable. This is possible since the aeroelastic interaction with the wake of the blade can give rise to complicated flow structures and also depends on the phase of the aeroelastic mode and the direction of the tower vortex. In Figure 9, we choose three different rotation instances to show the variation in the aeroelastic effect. It can be observed that aeroelasticity introduces a slight delay in the onset of the aerodynamic moment peak because of the tower wake.
Also, the peaks and range of aerodynamic moment are reduced especially in Figure 9A, while the smaller peaks introduced after the blade passes the tower is observed to be higher. In general, it can be said that blade flexibility results in smoother variation in the aerodynamic moment, which is more physically intuitive. However, the aeroelastic effect for the investigated cases in not very large, as the effect of the structural oscillations on the blade forcing is not significant in comparison to the effect of the tower wake. The aeroelastic investigations provided many interesting insights into the behaviour of downwind wind turbines. The aeroelastic code is able to reconstruct complex flow features, and we have demonstrated the ability of the code to match experimental measurements reasonably well.
However, because of the high computational costs, this technique is impractical for performing UQ. In the next section, we employ the code in a data-driven framework to develop ROMs in order to perform UQ.

REDUCED ORDER MODELLING FRAMEWORK
In order to perform UQ, we need inexpensive simulation codes. In this section, we introduce a framework for construction of ROMs for predicting aeroelastic forces. The framework is based on an earlier study by the present authors 37 applied to test case of the Goland wing; here, we discuss it concisely. It is based on building the ROM for prediciting the aerodynamics only, which is responsible for bulk of the computational time.
This ROM is then coupled to the structural solver to predict aeroelastic characteristics. The ROM is initially trained through aerodynamic data obtained from precribed forced motion to the structure. Two types of ROMs were developed in Sarma and Dwight 37 : an autoregressive with exogenous (ARX) model and a linear parameter varying (LPV)-ARX model. The ARX model is a time-marching method that predicts the output of a system based on earlier observations and inputs, represented in the form of a recurrence relation. An ARX model of a multi-input multiple-output (MIMO) system for output quantities Q j and exogenous inputs u 1 , u 2 , … , u N u at time t ∈ N can be represented by the following: ∀j ∈ {1, … , N Q }, where a and b are model coefficients to be trained. The lag operator time-shifts the quantities, q (t) ≡ h q (t−h) . Concisely in operator form, the system can be written as follows: where Q ∈ R N Q and u ∈ R N u are the N Q outputs and N u exogenous inputs, respectively. The operators  Q ( ) and  u ( ) are used to define the recurrence relation; the corresponding coefficient tensors A Q and A u can be written in the form: The elements of the tensor are estimated from a training phase, where time-domain data obtained from a high-fidelity solver (such as the one defined in the previous section) are employed. The model orders m and n define the number of past outputs and inputs considered. If we assume that the outputs are mutually independent, then many terms in A Q are zero, in particular: In the case of the wind turbine, u denotes the displacement of the blades in the fluid mesh, which is interpolated from f n in (5). The outputs Q j denote the pressure field on the surface of the blade. Furthermore, the ARX model is modified in order to account for the change in the blade moment when the blade is subjected to the tower wake, which introduces nonlinearity in the response. We introduce a localised forcing term G (t) in (8) to form the model: The localised forcing is in a Gaussian shape, which consists of a location parameter g, which locates the tower wake position, given by the following: where g defines the azimuthal position of the blade and is the width of the forcing. The parameter g can be used to parametrise the rotation speed of the turbine, which will be later incorporated in the UQ analysis. For training the model, the blade is prescribed a motion given by u †(t) ∈ R N u (exogenous inputs), and corresponding aerodynamic forces Q †(t) ∈ R N Q (outputs) are obtained from the flow solver, where † terms represent the training data. As usual, the training signal u † should be informative and plentiful in terms of frequencies and amplitudes of interest in order to identify A Q , A u , and A G . One obvious choice is to employ chirp signals for training; for the ARX model, a representative chirp is given by the following: , is known as the chirp rate with 1 , 0 , and Δt as the final frequency, initial frequency, and time-step, respectively. The prescribed amplitude 0 is in the order of millimetres for each excited mode. Independent modal excitations for each of the four modes are performed to obtain corresponding generalised aerodynamic forces.The excitation frequency for a particular system can be estimated based on the natural frequencies of vibration of the system, which can be obtained from the modal decomposition of the structure explained in Table 1. These natural frequencies are employed to tune ( ) in order to ensure that the frequency spectrum of the chirp overlaps with the natural frequency of the particular mode. As shown in Table 1, four structural modes are considered in this investigation. A statistical model is constructed given by the following: whereQ (t) ∈ R N Q is the predicted output at time instant t. The Gaussian error term (t) ∈ R N Q accounts for the discrepancy in prediction. A least-squares problem can be formulated to minimise (t) , given by the following:

Verification of ARX: Reconstruction of forced motion
Verification of the ARX model for reconstruction of forced motion is performed, where sinusoidal test signals are used to specify the structural motion of the blades. In the coupled aeroelastic solver, the modal force, also known as the generalised force Q n (t) shown in (5), is used. In Figure 10, the normalised generalised forces obtained from the full CFD solver and the ARX model are compared with TSR ≃ 1.6 and 6.3 of the wind turbine. It can be observed that the ARX solver is able to reconstruct the forces accurately close to the CFD predictions. The localised forcing term is also able to generate the trend in the peak encountered because of the tower wake. As explained in Figure 8, the change in TSR causes a change in the frequency and amplitude modulation of the low-frequency unsteadiness. It is seen that the localised forcing is able to reconstruct this behaviour with some discrepancy in the absolute value of the peaks; however, the trend is clearly estimated.

Verfication of ARX for aeroelastics: Estimation of blade forcing
After the verification of ARX through reconstruction of forced motion, we now employ the model to reconstruct unsteady aeroelastic forces. At fixed values of TSR and inflow conditions corresponding to experimental conditions, multiple ARX models are generated. The ARX models are coupled to the modal structural solver and aeroelastic loads on the blades are predicted. Figure 11 shows the estimation of the blade forcing at the two values of TSR. They are compared with CFD-based aeroelastic results. It can be seen that the ARX model is able to reconstruct the trend of the aeroelastic forcing accurately. For the lower TSR ≃ 1.6, we observe an offset in the absolute value of the mean forcing after the blade passes in front of the tower. However, at TSR ≃ 6.3, the offset is comparatively lower. This may be due to the introduction of the relatively stronger vortices from the tower wake ( Figure 6B) at the reduced TSR, resulting in higher mean forcing in the CFD solver. This effect is not captured by the ARX model, since the training phase involves forced motion of the blades, which is hence unable to capture this effect as it occurs only in the aeroelastic case. For an absolute agreement, the training of ARX has to be tailored by taking into account the phase shift between the forcing signal and the vortex generated by the tower. However, the model is able to estimate the forcing introduced by the tower accurately. In the next section, we use this ROM in order to propagate uncertainties in the experimental system.

UNCERTAINTY REDUCTION OF EXPERIMENTAL SYSTEM
The uncertainties existing in the investigated experiment are significant, and this can be seen in Figure 4B, where the structural moments obtained from the experimental measurements are shown at TSR ≃ 1.55. This experimental uncertain response is used as a reference for comparison to the uncertainty because of the structure s and rotation g, which will be introduced in this section. It should be mentioned here that in the experiment, 38 the distribution of the blade properties such as mass and damping were imprecisely known due to the material composition of the blade. The blades were manufactured with glass fibre and a foam core in the internal structure and surrounded by foam for the aerodynamic shape. The blade properties were estimated through an optimisation procedure by comparing the experimental and simulated eigenfrequencies.
In order to estimate these parameters in a probabilistic sense and define them in the Bayesian set-up, we introduce uncertainties in the parameters of (1) to (3), which gives an uncertain U n and n in (4) and (5), respectively. The structural parameters of the wind turbine that are considered uncertain are the distribution of mass, Young modulus and moments of inertia in two directions. Furthermore, as mentioned in Section 3, it is observed that the measurements of the rotational speed of the turbine are highly uncertain. Hence, the rotational parameter is also considered uncertain in this investigation, which is introduced through the term G in (10). As described in Section 4, G is defined by a location parameter g, which can be tuned to change the location of the tower wake. The components of g correspond to the rotational speed of the turbine for each rotation, which is considered to be uncertain here. The framework for the uncertainty propagation and reduction problem is based on an earlier investigation by the present authors, 37 which is briefly discussed below.

Methodology for uncertainty reduction
The effect of epistemic uncertainty in the input parameters s and g of a wind turbine on the generated structural moment is investigated. We use Bayes theorem 58 : where M is the experimental structural moment measurement. P 0 (s) and P 0 (g) are the prior distributions. P (M|s, g) is the likelihood derived from a statistical model relating the model prediction to the experimental observations, and it represents the difference between the experimental data and simulation model. The experimental noise is accounted here through a covariance matrix, generated from the experimental data. Finally, P(s, g|M) is the posterior, which updates the estimates on the input parameters s and g. The posteriors are then propagated through the simulation code to obtain a posterior predictive distribution of the structural moment.
Because of uncertain structural parameters, the structural modal matrix changes for each propagated sample. However, the trained ARX model (10) is based on mode shapes U s given in (4), obtained from the deterministic structural parameters s. Each sampleŝ drawn from the probability distribution will form a new mode shape Û s . As a result, both the generalised displacements f and the generalised load vector Q will change. For a new sample, the generalised displacementf has a new mode shape Û s . After transformation of displacements to the aerodynamic mesh, the corresponding modal deformation f in terms of U s to be used in the ARX model is solved through the minimisation problem Û sfs = U s f s . Similarly, the generalised load vectorQ forŝ to be used in the structural equation is obtained from where U a is the modal matrix for the aerodynamic mesh. This allows the trained ARX model based on the baseline model to be used for estimating the aerodynamic forces in the probabilistic framework.

Propagation of priors on the input parameters
The uncertainty on the inputs s and g is defined in order to obtain a probabilistic estimate on the structural moment. The cheap ARX model is used to obtain bounds on the structural moment by propagating these uncertainties. Lets andḡ be the mean (deterministic/nominal) values of the uncertain parameters. For a ±2% variation abouts and ±1% variation aboutḡ, we assume uniform probability distributions such that, s ∼  (0.98s, 1.02s), g ∼  (0.99ḡ, 1.01ḡ).
We assume a uniform prior, but the prior could be based on expert opinion. Also, the covariance matrices for the priors are diagonal, meaning we neglect correlations between the parameters. For each sample from the prior distribution, the numerical approximation of structural moment, which is a function of the load predictions from (10), can be written in the form: This provides us probabilistic estimates on the structural moment history. The Bayesian framework is now introduced, which is utilised to identify the input parameters from data.

Bayesian framework for uncertainty reduction
Let M + ∈ R n be a measured value of the structural moment vector; n refers to the azimuthal positions at which measurements are obtained.
We now introduce an additive random variable to account for noise in the measurements, specified as an unbiased normal distribution with a covariance matrix , obtained from the experimental data (variance due to the moments reported for repeated runs or rotations of the turbine).
Neglecting systematic modelling errors as in Sarma and Dwight, 37 the statistical model can be written as follows: which gives the probability of observing M + given s and g, also known as the likelihood as follows: P(M + |s, g) ∶= P (M + −M(s, g)).

Expression for posterior of input parameters
The priors on the structural parameters can be now updated with the likelihood given by (18) using Bayes theorem to obtain the posterior 58 : P(s, g|M + ) ∝ P(M + |s, g)P 0 (s)P 0 (g), ] P 0 (s)P 0 (g). The standard checks to obtain convergence in MCMC were applied. 59,60

Propagation of input uncertainties
The use of the ARX model enables computationally cheap aeroelastic estimations; thus, the uncertainties can be propagated with brute force Monte Carlo method. We first investigate the effect of uncertain s on the structural moment of the turbine. Thereafter, the effect of uncertain g is ascertained, and finally, all the uncertainties are combined to form a realistic prior for the Bayesian identification problem. Figure 12 shows the confidence intervals of the structural moment at 30% span location. The confidence intervals are plotted for 1%, 2%, and 5% variation (shown by the gray shaded area) in the structural parameters with respect to the mean (deterministic) values (shown by the black solid line).
The variation in the structural moment is introduced by the change in the aeroelastic frequency of vibration of the blade because of uncertain s.
This effect is however neutralised as the blade approaches the tower and is subjected to the wake forcing. This can be clearly observed in responses from higher variability in the input forcing, such as Figure 12B. The tower forcing is introduced at around 1.5 azimuthal position. The uncertain response has a staggered distribution around the mean, which can be observed from the crests and troughs of the harmonic response.
The shaded area (uncertain response) in the crest after the 1.5 azimuthal position is symmetric about the mean in terms of the frequency, while FIGURE 12 Confidence intervals of structural moment at 30%-span for TSR≃1.55 assuming uniformly distributed input parameters: blade sectional distribution of mass, Young modulus, and moments of inertia I and I. TSR, Tip Speed Ratio FIGURE 13 Confidence intervals for structural moment at 30%-span location at TSR≃1.55 assuming uniformly distributed localised forcing term in the ARX model. ARX, autoregressive with exogenous; TSR, Tip Speed Ratio the variation is only in terms of difference in amplitude. The effect of the tower wake thus collapses the frequency change introduced by the uncertain inputs. For the case with 5% variation (Figure 12C), the uncertain responses for structural moment have a significant spread because of the larger change in structural frequency. This large uncertain response with 5% variation is considered to be large for the current investigation.
In the next step, rotational uncertainty is introduced by considering g to be uncertain. Figure 13A,B shows the uncertain response in the structural moment for 1% and 2% variation in g. It is observed that the uncertainty in the location of the forcing term introduces larger spread in the response of the structural moment after the 1.5 azimuthal position. This is expected since the localised forcing term takes into account the forcing from the tower wake and the rotational uncertainty would in turn result in a phase shift in the incidence of the tower wake. When 2% variation is considered ( Figure 13B), the spread is comparatively large leading to significant difference with respect to the experimentally observed uncertainties.
In the final propagation problem, both the structural and localised forcing uncertainties are considered; the corresponding structural moment response is shown in Figure 14. As already discussed, for the localised forcing term, 1% variation is considered. For the structural parameters, both 1% and 2% variations are considered. The structural moment response is found to be symmetrical about the mean location. In order to consider for the higher spread in the response around the 1.5 azimuthal position, we consider 2% and 1% variation in s and g respectively in the Bayesian identification problem.

Reduction of uncertainties
Experimental data (shown in Figure 4B) is now introduced in (19) in order to identify the input parameters of the wind turbine. As already discussed, in this section, we consider 2% and 1% variation in s and g respectively. Figure 15 shows the comparison of the considered uncertainty with respect to the mean of the experimental measurements. In the Bayesian reduction problem, there are two considerations in the analysis presented here. First, the aeroelastic uncertain response has significant overlaps with respect to the experimental mean in most azimuthal positions, except near the 1.5 azimuthal position. Around this point when the blade is in front of the tower, the discrepancy is large. Second, there is large variance in the experimental measurements, as shown in Figure 4B. Based on this, posteriors are identified considering full and FIGURE 14 Confidence intervals for structural moment at 30% span location at TSR ≃1.55 assuming uniformly distributed localised forcing term and structural parameters: mass, Young modulus, and moments of inertia I x and I y in the ARX model. ARX, autoregressive with exogenous; TSR, Tip Speed Ratio FIGURE 15 Experimental mean and structural moment response for 1% and 2% variations in localised forcing and structural parameters FIGURE 16 Posterior distributions with 1-D histograms of structural parameters s 1 and s 2 (mass and Young modulus), localised forcing parameters g 1 to g 4 and 2-D scatter plots for each parameter pair truncated azimuthal position history. For the truncated analysis, the data around the 1.5 azimuthal location is ignored. Furthermore, the experimental variance is reduced to ascertain the effect on the identification.
When the full azimuthal position history is considered, the structural parameters are not identified, while localised forcing parameters are somewhat identified. However, the reduction does not reveal a clear correlation, even with reduced variance, which is attributed to the large difference between the experimental and aeroelastic estimates around the tower forcing location. Here, we would like to refer to our earlier discussion in Section 3.1 regarding the source of uncertainty in the tower location. This could be attributed to the limitations of the aeroelastic model. This discrepancy is avoided in Figure 16, where the posteriors identified with the truncated azimuthal position data are shown. We observe that in this case, both s 1 and s 2 are identified, while g 3 and g 4 are also skewed in reference to their uniform priors and they show a clear negative correlation in the 2-D histogram. Similarly, s 1 and g 4 are also scattered towards the lower edge of the prior uncertain domain. In the next step, we propagate the posterior distributions to predict the posterior predictive distribution of the structural moment. Figure 17A shows the confidence intervals under prior distribution for the assumed uncertainty, while the confidence intervals with identified posteriors (Figure 16) is shown in Figure 17B, which clearly reveals a reduction in uncertainty at all azimuthal locations. In order to obtain a quantitative comparison, we select time or azimuthal position instances (shown in Figure 17A) and plot the probability density estimates at these locations.
The time instances are marked t 1 to t 6 in Figure 17A, which are selected such that structural moments at crests, troughs, and between them are obtained, with more time locations considered around the 1.5 azimuthal position. From Figure 18, it is observed that there is reduction in the variance with respect to the prior at all locations. It can be seen that there is a larger reduction in uncertainty at locations t 2 and t 4 , compared with other instances. Both these time instances correspond to azimuthal positions located between the crest and trough, at which the identification data was used in the likelihood. At t 6 , which is after the blade passes the tower, the reduction in variance is comparatively less, and a large uncertainty is still observed. This is expected since the experimental measurements around this location were not used to update the posterior, and hence, it is not informed.  Probability density estimates of structural moment at 30% span location at TSR≃1.55 for t 1 -t 6 . TSR, Tip Speed Ratio To summarise, it can be said that because of the large uncertainty in the experimental data and the inadequacy of the model to provide accurate predictions near the tower location, consideration of data near these locations are unable to reduce the parametric uncertainties significantly.
The uncertainty because of the forcing from the tower wake can lead to large uncertainties in the structural moment, and it is attributed to the modelling uncertainty here, primarily to the RANS-based fluid solver, which may not be able to predict blade loads accurately in comparison to LES or DNS. As mentioned earlier in Section 2.2, the computational requirement of the latter solvers is not feasible for the ROM training.
Because of the complex nature of the wake dynamics, the discrepancy between the experimental uncertainties and aeroelastic responses is large at these locations. As a result, the reduction in uncertainties of the structural moment after identification is not large. Despite this, uncertainties have reduced at all azimuthal positions as evident from Figure 17. Also, it has been established that the localised forcing terms, which take into account the tower wake, can be identified with this methodology (Figure 16). This uncertainty can be critical in terms of evaluating the fatigue behaviour of a wind turbine blade under uncertain operational conditions. For future developments, model uncertainty can be introduced in order to account for the discrepancy near the tower.
In this section, the effects of the considered structural and rotational uncertainties are shown. The impact of these uncertainties are significant in the operation of the turbine. The structural uncertainty affects the vibration characteristics of the blade and thus the blade dynamics. This also could affect the stability characteristics and fatigue of the blade. A probabilistic estimate of these quantities of the blade can be useful to reduce risks during operation of the turbine. The other uncertainty considered is the rotation speed of the blade, which is a likely scenario in an operational wind farm because of changes in wind speed, direction, and the tolerances associated with mechanical components such as gear.
This has a direct impact on the power output of the turbine. Again, a probabilistic estimate on this parameter can assist wind farm designers to estimate power output accurately.

CONCLUSION
In this paper, the aeroelastic behaviour of an experimental wind turbine system is investigated with a full detailed 3-D aeroelastic model considering all structural components. The developed aeroelastic model is shown to be able to provide predictions within the bounds of the experimental structural moment measurements. A ROM is developed to predict aeroelastic characteristics using ideas of system identification and a localised forcing term, which enables computationally cheap estimation of structural moments. It is shown that the model can provide accurate structural moment history in comparison with CFD-based full aeroelastic solver. Furthermore, uncertainties in structure and rotational speed of the turbine are considered in order to provide probabilistic estimates on the aeroelastic behaviour. The uncertainties are reduced using a Bayesian framework and incorporating the experimental data in the likelihood estimation. It is established that the methodology can be used to reduce uncertainties in aeroelastic estimates. The potential applications are predicting fatigue in wind turbine blades, which can be important to consider for downwind oriented wind turbines because of the incidence of tower wake. The work can be extended to characterise the low-frequency unsteadiness of the tower wake, which can be included in the formulation of the ROM. Furthermore, better aeroelastic estimation of blade forces near the tower with respect to experimental measurements can improve the identification of the input parameters of the wind turbine.