Abstract
Rapid development of offshore wind foundation models has resulted in a large number of built structures with generally underestimated foundation stiffness properties and a need to update and validate both the individual structural models and the underlying foundation design frameworks. This paper outlines a structural health monitoring approach, based on the combination of output only structural health monitoring methods and model updating, to estimate foundation stiffness parameters using field monitored data. Field monitoring data from an offshore wind turbine under idling conditions, over a large monitoring period, are presented and operational modal analysis is applied to estimate the modal parameters. Those are compared to modal properties predicted by finite element models, employing either old (API/DNVGL) or new (PISA) foundation design properties, which are calibrated using geotechnical site investigation data. A new approach to interpret seabed level statically equivalent foundation stiffness, in terms of effective lateral and rotational stiffness against load eccentricity, is presented. Seabed level statically equivalent foundation properties are updated by comparison against the observed modal behaviour and the optimised foundation parameters are presented, demonstrating a close match to the predictions of the PISA method.
1. Introduction
There is currently over 40 GW globally of installed and ageing offshore wind turbines (OWTs). Almost all of these structures were designed using now obsolete methods, such as the traditional API/DNVGL monopile foundation p–y design approach. Monitoring of offshore wind structural natural frequencies [1] and recent foundation model development (e.g., PISA) suggests significantly higher stiffness and capacity than previous design standard models. There is a pressing need to reassess these structures, for asset valuation, remaining life assessment, and damage identification [2] and to support further development of design models.
Already widely used in drive train condition monitoring, structural health monitoring (SHM) approaches are being developed to inspect and assess the wider range of superstructure components. These approaches span from input-state-parameter estimation methods, which enable estimation of both the properties and input excitations [3], through to operational modal analysis (OMA) methods [4], which have been used to infer damage or change from deviations in structural modal properties[5]. The application of a number of OMA methods to OWT monitored data has been explored by different groups, such as the work of Oliveira et al. [6], which showed that stochastic subspace identification (SSI) methods can provide robust estimates of OWT modal properties over a range of operational modes. Such methods can be extended to directly identify structural properties using a second stage of model updating, where numerical structural model parameters are optimised to best match the predicted and measured modal behaviour. Whilst model updating has been used successfully to identify foundation properties in applications such as bridge, breakwater, and axial pile analyses [7–9], challenges in validating measurements, due to a lack of accurate modelling or benchmark data, have prevented implementation in directly estimating OWT foundation properties from monitored data. Whilst Anderson [10] investigated the components of a framework for OWT foundation parameter estimation, this was explored only through numerical simulation and applied to synthetically generated data.
Poor prediction of short rigid OWT monopile foundation behaviour, using existing methods developed for long slender oil and gas pile geometries [11], led to a number of research programs to develop new offshore wind monopile design methods [12–15], which promise significant reductions in design conservatism and cost. These methods are typically based on numerical foundation modelling or scaling of laboratory element testing and while some have been compared against reduced scale laboratory or field testing [16–18], they have not been validated against full scale structures. The lack of field validation introduces the risk of underconservatism and increased sensitivity to effects, such as cyclic degradation or scouring.
This paper describes a model updating SHM approach for the estimation of OWT monopile foundation model parameters, by calibrating the modal properties of a finite element model to those estimated from data. Section 2 presents the monitored field data, followed by the application of an OMA approach to accelerometer measurements and then the identification of stabilised modes for a full scale OWT monitored by Parkwind and VUB. Section 3 outlines the adopted numerical models, including a description of the finite element formulation, the parameters for both a flexible and rigid-rotor model, and the models adopted for predicting the soil-foundation macrostiffness. Section 3 also derives and demonstrates a novel approach to present seabed level statically equivalent foundation stiffness in terms of effective lateral and rotational stiffness against load eccentricity. Section 4 describes and applies a model updating approach, to update the properties of the statically equivalent seabed level foundation macro-element and compares the resulting optimised foundation behaviour against previous (API/DNVGL) and new (PISA) design methods.
2. Field Data
2.1. Wind Turbine Description
Data come from an instrumented Vestas V112-3.3 MW offshore wind turbine, in the Nobelwind farm located 47 km from the coast of Belgium in the North Sea. The geometry and accelerometer sensor layout are illustrated in Figure 1. The structure is mounted on a 5.0 m diameter monopile embedded 29.9 m into a layered stratigraphy of dense sand overlaying overconsolidated clay. The mean lowest astronomical tide (mLAT) is 31.7 m above seabed level. A transition piece connects the monopile through a grouted sleeve connection and then to the tower through a bolted flange. The turbine has a diameter of 112.0 m with a hub height of 110.5 m above seabed level.
2.2. Monitoring System
Perpendicular biaxial accelerometers measure at three tower levels: below the transition piece-tower flange (ACC-1); the midtower service platform (ACC-2); and below the main yaw bearing (ACC-3). Accelerometer signals are sampled at 3200 Hz and downsampled to 100 Hz for logging of continuous 10-minute datasets. The nacelle mounted supervisory control and data acquisition (SCADA) system logs 10-minute average data, including wind speed at the rotor axis elevation, blade pitch, rotor speed, and yaw angle. Raw tower accelerometer signals are transformed to fore-aft (FA) and side-side (SS) components using the 10-minute averaged SCADA rotor yaw angle. Wave radar measurements, 2.5 kilometres from the turbine, provide 10-minute average measurements of tidal elevation.
The model updating approach in this paper focuses on datasets that come from two periods of low-wind idling turbine conditions, indicated with grey shading in Figure 2, totalling 87 10-minute datasets. The idling periods feature no change in yaw angle across consecutive datasets, so that the risk of error in the acceleration yaw transformation is mitigated.
(a)
(b)
(c)
2.3. Stochastic Subspace Identification
In the absence of measured load excitation, the subspace state-space system identification method N4SID [19] is applied to the fore-aft and side-side accelerometer signals to identify the modal parameters of the system from the monitored data.
The full algorithm is presented in Van Overschee and De Moor [19] and a brief description of the method follows. A discrete state-space system is described in terms of the state, , input, , output, , and feedthrough, , matrices describing the transition of the system between steps and , in state space form as follows:where and are the process noise and measurement noise vectors, respectively. The steps correspond to a finite sampling time step . The sampling frequency used after downsampling the data was 100 Hz. The method requires the discrete input, , and output, , data to identify all state-matrices, or alone to identify and , as is the case in this study. is then converted to the state matrix of the equivalent continuous state-space system usingwhere is the matrix logarithm.
The eigenvectors and eigenvalues λ are obtained by solving the eigenvalue problem:
The matrix is used to transform the eigenvectors from the state vector, z, basis to the output vector, y, basis (sensor basis):
Mode shapes presented from here on are in the sensor basis and the presuperscript s is omitted for brevity.
25 seconds are removed from beginning and end of each dataset to reject any filtering-related distortion.
The procedure described in equations (1)–(4) is repeated for different model orders, to obtain the eigenvalues of the corresponding for that model order. A stabilization procedure is then followed [4]. A pole is considered stable if it presents a variation of frequency smaller than 1% and a modal assurance criteria (MAC) parameter mode shape similarity of 0.97 or greater relative, to a pole from the previous model order. Poles with damping ratio of 20% or more are rejected. The frequencies corresponding to the poles identified by the N4SID algorithm are illustrated in Figure 3 for a representative dataset and model orders up to 300. Identification to high model orders is necessary to capture the influence of unmeasured excitations and to ensure stabilisation of the modes of interest [4]. Frequencies presented in this paper are normalised, , by the identified frequency of the 1st tower dominated side-side mode, as described in Section 2.4, averaged across all datasets. Modes are classified “stable” if they are both frequency- and MAC-stable.
2.4. Tower Dominated Mode Shape Identification
The model updating algorithm described in this paper utilises an objective function that is based on the difference between the numerically predicted and the N4SID estimated tower dominated mode frequencies. The following procedure is used to determine the tower dominated modes for a given dataset. If the system were proportionally damped, we might expect to identify a pair of orthogonal lateral modes comprising the nth pair of tower dominated frequencies, one fore-aft and one side-side. In practice, measured mode shapes exhibit motion due to contributions from multibody dynamic coupling and nonproportional damping, both predominantly from the blades, which means that the real part of the identified mode shapes is unlikely to align with the principal structural axes. Numerical simulations by Anderson [10] demonstrate that nonproportional damping occurs due to the aeroelastic interaction of idling rotor blades, resulting in significant deformation of the observed mode shape from the undamped mode shape. In agreement with the work by Sinha [20], Anderson showed that the mode shape influence from nonproportional damping can be mitigated, by approximating the mode shape corresponding to the underlying undamped system, using a “complex normalisation” technique, as originally presented by Fillod [21].where j is the element of the complex mode shape that maximises .
“Complex normalisation” is in contrast to the conventional “real normalisation,” which assumes predominantly proportional damping, allowing the imaginary component to be neglected, i.e.,
Figure 4 presents the mode shapes of the 1st, 2nd, and 3rd pairs of tower dominated modes in the fore-aft and side-side directions, with both “real” and “complex” normalisation, for a representative dataset. Close spacing of the 1st pair of tower modal frequencies prevents delineated identification of mode shapes that align with the principle fore-aft and side-side axes and the mode shape direction shows significant variation throughout the idling datasets. However, comparison using the S2MAC [22] of the two identified 1st tower mode shapes against the corresponding numerical closely spaced modes, predicted in Section 4, indicates that the identified modes fall within the subspace spanned by the numerically predicted mode shapes.
Complex normalisation has little influence on the mode shapes of the closely spaced 1st tower dominated modes. However, the significant nonorthogonal mode shape components in the 2nd pair of real normalised tower modes are reduced when using the complex normalisation. This reduction was consistent across the available idling datasets and agrees with similar observations by Anderson [10] from simulated data. The complex normalised mode shapes show the expected characteristic shape and alignment towards principle axes. The side-side mode shape occurs at a higher frequency, due to the lower mobilised in-plane rotational inertia of the idling rotor, when compared to the fore-aft plane. 3rd tower modes exhibit the expected classical mode shape, but only the fore-aft mode shows differences between the real and complex mode shape normalisation.
A semiautomated approach is adopted to assign identified modes with tower dominated mode classifications. Manual inspection of stable modes is used to determine bounds of expected tower dominated mode frequency and an initial empirical estimate for each mode classification, as shown in Table 1. For modes that fall within these bounds, the dominant direction of the complex normalised mode shape is used to classify as either a fore-aft or side-side mode. The 1st tower mode pair is an exception to this rule, as the closely spaced frequencies result in the identified mode shapes being a linear combination of the corresponding fore-aft and side-side mode shapes. In this case, they are classified based on frequency sorting, with the fore-aft being the lower of the two identified frequencies. Here, multiple suitable modes fall within frequency bands, the closest to the initial empirical frequency estimate (see Table 1) is selected.
Figure 5 shows the classification for the stable modes across the 87 idling datasets. The classified 1st tower fore-aft and side-side mode frequencies are intermingled, which demonstrates the source of the difficulty for the N4SID algorithm to delineate the two orthogonal modes. The 2nd tower mode frequencies vary across the datasets, due to variability of water added mass throughout the tidal cycle. Higher sensitivity of 2nd tower modes to effects close to the seabed, such as scour or tides, has been well documented [5, 23]. Figure 5(c) also plots the 2nd tower frequency predicted by the optimised numerical model outlined in Section 4, which demonstrates good predictability of the frequency variability with the tide over the two continuous sets of idling data. Figure 5(d) indicates that the 3rd tower frequencies are relatively insensitive to changes in tidal level.
(a)
(b)
(c)
(d)
Table 2 presents the mean and standard deviation of damping estimates for the identified tower dominated modes across the 87 idling datasets, showing consistent values and variance to similar OWT SHM campaigns [24].
3. Numerical Structural Modelling
A numerical model simulates the OWT dynamic response and the foundation soil-structure interaction and is used as the basis for the subsequent model updating foundation parameter estimation in Section 4.
Two foundation model variants are used, one with a finite element representation of the elastic monopile and depthwise-varying soil reactions and one with a macroelement to represent the overall elastic response of the foundation at seabed level. The depthwise-varying foundation model is used to provide design estimates of the expected foundation behaviour, whilst the macroelement model is later used as a basis for the foundation parameter estimation model updating, described in Section 4.
Two variants of the turbine superstructure are modelled, one where the flexibility of the blades is taken into account, referred to as the flexible-bladed (FB) assembly and one where the blades are assumed to be rigid, named the rigid-bladed (RB) rotor assembly. The flexible assembly is highly informative in developing an understanding of the system dynamics and choosing the model updating parameters. However, it is common for detail of the rotor-nacelle-assembly structure to be unavailable to substructure designers, and so, a model updating approach is sought based on the simpler-rigid-assembly model. Both models and their predicted behaviour are described in the following sections.
3.1. Finite Element Modelling
The wind turbine was modelled using a Timoshenko beam finite element model. Modelling of shear compliance is necessary for accurate prediction of embedded pile stiffness [25] but is not anticipated to play a significant role in the dynamic modelling of the slender tower. The formulation of the lateral components of the stiffness and mass matrices follows the approach outlined by Friedman and Kosmatka [26], extended to biaxial lateral deformation. Mass and stiffness matrices are evaluated over each element using Gauss integration with four gauss points per element. The Friedman and Kosmatka approach derives transformed element shape functions that exactly satisfy the homogenous form of the static equations of a uniform Timoshenko beam. Control of a single parameter, β, the ratio of the beam bending stiffness to the shear stiffness as described in equation (8), incorporates the influence of shear compliance and can be reduced to the stiffness and mass matrices associated with conventional Euler–Bernoulli beams for .
Conventional cubic Hermite interpolation functions are adopted as the basis of the lateral degree of freedom shape function formulation, in a similar fashion to Friedman and Kosmatka [26], and linear Lagrangian interpolation functions are adopted for axial and torsional degrees of freedom [27]. The resulting structural and soil stiffness elements provide identical lateral behaviour to those described in Byrne et al.’s study [25].
12 elements model the embedded pile and 49 elements model the monopile stick-up, transition piece, and tower design geometry. For the flexible rotor model, blades are simulated with 10 Euler–Bernoulli elements per blade. Monopile, tower, and nacelle properties are modelled using characteristic design geometry and parameters, with the properties of the transition piece beam elements derived through 3D finite element modelling. It is acknowledged that inaccuracy in the as-built design will influence the subsequently identified foundation properties, but it is expected that such errors are small relative to the uncertainty in the modelled soil-structure interaction stiffness. All steel components are modelled with Esteel = 205 GPa. The density of the monopile is 7850 kg/m3, whilst the density of the transition piece and tower are 8500 kg/m3, to account for paint and ancillary steel. The grout sleeve between the monopile and transition piece is modelled with Egrout = 40 GPa and a density of 2400 kg/m3. The distributed added mass per unit length allocated for the internal and external water is given by the following equations:where Cm is the added mass coefficient of the external fluid. For this study, Cm = 0.5 was found to best match the range of 2nd tower natural frequency variation with tidal level over the observed period. The level of the internal water is assumed to coincide with the external free surface.
3.2. Rotor-Nacelle-Assembly Modelling
Available information on the rotor-nacelle-assembly structural properties of the monitored OWT are commercially sensitive and limited to the mass, second moments of inertia, and centre of gravity offset from the tower top, which provides sufficient data to formulate only a rigid-rotor-nacelle-assembly model. However, additional insight into the system dynamics can be gained through simulation of a flexible rotor-nacelle-assembly.
A flexible rotor model has been approximated through scaling of the benchmark NREL 5 MW wind turbine rotor [28]. Blade properties and geometry are linearly scaled based on the ratio, , of the diameters of the NREL 5 MW (126 m diameter) and the Vestas V112-3.3 (112 m diameter) rotors. A subsequent scaling factor is applied to the blade mass per unit length, , to match the rotor moment of inertia to that of the V112-3.3. This further reduction is expected, considering the lower rotor design power take-off density, 335 W/m2, for the Vestas V112-3.3 turbine, compared to the 401 W/m2 of the NREL 5 MW, by a ratio of . Lower power density will result in lower thrust density and a lower blade force per unit length. The blade root moment, which is a common design driver in blade design, therefore scales with . To preserve similarity of blade root stress, the blade section second moment of area and moment of inertia are also scaled by . The overall scaling of each blade structural property is the product of the geometric, mass, and power take-off scaling factors, as shown in Table 3.
In order to predict the modal properties of the turbine with an idling rotor, the rotor shaft is modelled with negligible torsional stiffness and relatively rigid flexural stiffness.
3.3. Elastic Foundation Modelling
It is assumed that the foundation will see small deformations during idling rotor conditions, which can appropriately be modelled using linear-elastic foundation models. Foundation soil-structure interaction is modelled using two approaches. The first approach is shown in Figure 6(a), which simulates the embedded monopile and applies depthwise-varying linear soil reactions and concentrated base reactions axi-symmetrically in the x- and y-directions, using the method outlined in Byrne et al.’s study [25]. The second approach models the macro embedded foundation response through axi-symmetric seabed level coupled linear springs, KG, as shown in Figure 6(b) and equation (11), where subscript G denotes ground level. The influence of the embedded foundation mass and inertia is assumed to be negligible. Vertical displacements and torsion are constrained at the ground level.
(a)
(b)
3.4. Depthwise-Varying Foundation Modelling
Depthwise-varying soil reactions are calculated using soil properties derived from the wind farm consultant Geotechnical Investigation Report, which included geotechnical site investigation interpretation at the turbine location. Note that the site investigation and consultant design focussed on calibration of a typical API/DNVGL p–y foundation design model. Table 4 outlines the depths and characterisation of the soil units. Figure 7 plots the assumed profiles of soil properties. Scour protection is conservatively assumed to provide no lateral reaction and only provide overburden to the soil vertical effective stress. Detailed small strain shear modulus measurements in the shallow sand layers are unavailable, so the Hardin [29] stiffness profile is adopted, as shown in equation (12). Correlations between relative density and voids ratio, and the value of B = 875, matched the PISA project Dunkirk sand model [30], which is assumed to closely approximate the sand characteristics at the nearby Nobelwind farm.where G0 is the small strain shear modulus, B is an empirical constant, p’ is the mean effective stress, e is the voids ratio, and pref is a reference pressure of 101.3 kPa. G0 in the clay layer (below 8 m depth) is approximated from an average profile of highly scattered seismic cone penetrometer measurements from a nearby borehole. The profile of undrained shear strength, su, (Figure 7(c)) is based directly on the consultancy report values.
(a)
(b)
(c)
Depthwise-varying soil reactions are calculated using either the PISA method, outlined by Byrne et al. [25] and Burd et al. [31] for clay and sand layers, respectively, or the API-DNVGLp-y method [32]. In this study, the implementation of the API-DNVGL approach for sand is modified to calculate the depth below seabed based on the equivalent depth, calculated based on the local unit weight and vertical effective stress, to account for the influence of scour protection overburden:where C1, C2, and C3 are coefficients dependent on the friction angle and
Note that z is defined as negative below ground here.
Kallehave [1] presented a modification to the sand subgrade modulus predicted by the API/DNVGL method, using equation (15), which is now well recognised within industry practice. The use of this modification is also explored and referred to as API/DNVGL-Kallehave.where k is the API/DNVGL predicted subgrade modulus, z0 = −2.5 m, D0 = 0.61 m, and m = 0.5.
3.5. Macroelement Foundation Modelling
Assuming predominantly elastic foundation response under idling rotor conditions, the depthwise-varying-foundation response at ground level can be modelled using statically equivalent coupled linear springs, called the macroelement. The statically equivalent macroelement is calculated from the depthwise-varying-foundation stiffness matrix Kdv, which is initially inverted. The ground level degrees of freedom are then extracted (subscript Gdof) and again inverted to provide the ground level macrostiffness matrix KG, as shown in the following equation:
The macroelement components for each of the PISA, API/DNVGL, and API/DNVGL-Kallehave foundation model types are shown in Table 5.
3.6. Equivalent Lateral kH and Rotational kM Stiffness
It is difficult to compare the behaviour of two coupled macroelements, as significant changes in individual elements of the stiffness matrix can compensate to provide minimal change in overall behaviour. The following analysis describes the derivation of the resultant lateral and rotational stiffness of the macrofoundation response and outlines an approach for graphical comparison of the elastic response. A transformed set of macroelement parameters are presented that are both an improved measure of foundation response and improved optimisation parameters for the subsequent model updating process, described in Section 4.
First, we define the ratio of shear force and moment at seabed level through a static equivalent vertical eccentricity h, of an applied lateral static load:
The equivalent lateral kH and rotational kM stiffnesses of the coupled linear spring can then be derived as a function of the load eccentricity:where and ψG are the ground level displacement and rotation.
Note that both kH and kM are functions of all three stiffness terms k11, k22, and k12, as well as the load eccentricity, h. Figure 8 illustrates typical variation of kH and kM against h. Three pertinent characteristics of the equivalent stiffnesses emerge, which are the intercept of the lateral stiffness, , the asymptote of the rotational stiffness, , and the negative gradient of the lateral compliance, :
(a)
(b)
(c)
These components can also be derived based on the components of the macroelement compliance matrix and are referred to as the compliance terms:
Note that whilst all the three macrostiffness terms, k11, k22, and k12, influence kH and kM, the compliance term C11 only influences kH and C22 only influences kM, as shown in equations (24) and (25). The individual terms of the compliance matrix therefore provide a more direct indication of the difference in foundation response between two macroelements.
Compliance term values for each of the three foundation models are presented in Table 5 and the corresponding variations of kH and kM against h are plotted in Figure 9. The basic API/DNVGL model and augmented subgrade modulus of the API/DNVGL-Kallehave model predict lower kH and kM than the PISA foundation model. This matches previous observations of natural frequencies in the field [1] and design model comparisons [33, 34].
As shown in Figure 9, kH and kM are highly dependent on the equivalent load eccentricity h. Whilst estimation of the load eccentricity applied by the quasistatic action of wind or waves is relatively straightforward, the equivalent load eccentricity during modal excitation should also be considered. The modal equivalent load eccentricity hi, corresponding to the mode specific applied moment MG,i over the horizontal modal force HG,i, may be estimated from a FE modal analysis. For any deformed mode shape of the structure the ratio of the ground level displacement, and rotation ψG,i,
Combining equations (11), (17), and (26) results in a matrix relationship between the applied loading and the ground level deformation, in terms of hi and ηi,
Dividing the second row of this matrix equation by the first results in an expression for the modal equivalent load eccentricity as a function of the stiffness terms and ηi,
The values of hi calculated from finite element modal analyses, featuring each of the design foundation models, are shown as circles in Figure 9. The load eccentricity associated with the mode shape of the 1st SS, 2nd SS, and 3rd SS tower dominated side-side mode shapes is plotted, but values are similar for the fore-aft mode shapes. Here, we note that hi for the dynamic behaviour of a given mode shape decreases as we move to higher tower modes.
3.7. Flexible and Rigid-Rotor Modelling Comparison
Figure 10 compares the modes identified by the N4SID algorithm for a representative dataset, against modes predicted by the rigid-bladed (RB) and flexible-bladed (FB) rotor variants of the finite element model. Figure 11 performs a similar comparison of the mode shapes of the 1st, 2nd, and 3rd pairs of tower dominated modes in the fore-aft and side-side directions. Foundation reactions are modelled using the PISA method, which is later shown in Section 4 to provide a good estimate of the optimal foundation stiffness.
Figure 10 demonstrates the identification of a large number of modes, some associated with the typical fore-aft and side-side tower dominated behaviour, as well as blade and rotor dominated modes. In practice, all identified modes feature contributions from several structural components through dynamic coupling. However, the rigid-bladed FE model has no degrees of freedom associated with blade flexure and therefore explicitly predicts only tower modes, which are a small subset of the stabilised modes shown in Figure 10. Similarly, the corresponding mode shapes shown in Figure 11 have no influence from multibody coupled dynamics and show perfect alignment with the principal axes. Table 6 lists the subset of modes from the identification and numerical modelling that have been identified as the tower dominated modes, which indicates that the rigid-bladed model generally achieves good predictions of the identified frequencies, except for the 3rd tower fore-aft mode. The rigid-bladed model also achieves good MAC comparisons to the higher identified modes, with low MAC values for the 1st tower modes, reflecting the variable direction of the mode shapes identified from the monitored data in this case, shown in Figures 11(a) and 11(c).
The flexible-bladed model predicts a greater number of frequencies than RB, which corresponds with data identified modes shown in Figure 10, such as the four identified modes around the 2nd tower dominated mode frequencies (). In some cases, predicted modes are not identified, such as some of the 1st blade flapwise dominated modes around ≈ 2.8. This is likely to occur due to tower accelerometers being insufficiently excited to make blade modes practically observable. Figure 11 shows that multibody dynamic coupling results in the mode shapes of the FB model not being perfectly aligned with the FA and SS planes. However, geometric and property differences between the scaled NREL benchmark turbine and real turbine result in out of plane components for the 1st and 2nd tower modes that differ from the observed data across the available idling datasets. Despite potential discrepancy of the rotor model, a good match is achieved between the majority of the predicted and identified tower dominated mode frequencies, with the exception of the 3rd tower fore-aft frequency as observed in Table 6. The inability to predict this particular mode with both rotor models is attributed to a particularly high sensitivity of the fore-aft tower-blade coupling at this frequency. For such a mode, accurate knowledge of the properties of the blades becomes more important. During idling, there is a full moment transfer between the rotor and tower in the fore-aft direction, as illustrated in Figure 12(a), but there is a reduced transfer of moment about the rotor axis and there will be a reduced influence of rotor-tower coupling effects in the side-side direction, as illustrated in Figure 12(b). This makes modes in the side-side direction better candidates for the model updating, given the uncertainty about the properties of the rotor-blade assembly.
(a)
(b)
To further support this observation, Figure 13 plots the variation of the natural frequencies predicted by the flexible rotor finite element model, featuring a PISA foundation, where the blade stiffness properties (E and G) have been multiplied by a scaling factor α. An α value of 1 indicates the expected blade flexibility and increasing α represents the transition to a rigid-rotor assembly at α = ∞. Figures 13(b)–13(d) show focussed plots of the 1st, 2nd, and 3rd tower dominated mode frequencies, respectively. As blade stiffness increases, there is an increase in the majority of the observed natural frequencies. However, some frequencies are relatively invariant with blade stiffness, which appears as pseudohorizontal bands, such as at ≈ 11.3, as shown in Figure 13(d). These α-insensitive frequencies (for α ≥ 1) include the 1st tower mode pair and those higher modes associated with the side-side tower dominated motion of the structure. The 1stfore-aft dominated tower mode is also relatively invariant with blade stiffness, as the 1st blade natural frequencies are significantly higher and there is negligible dynamic interaction.
(a)
(b)
(c)
(d)
The subsequent model updating will adopt the modes that show minimal frequency sensitivity to blade stiffness variation and which will minimise the risk of identification error due to uncertainty in the blade stiffness properties. This is particularly applicable to the common case where detailed blade design information is not made available by the turbine supplier to the substructure designers. These frequencies can be robustly predicted using a rigid RNA, where the prediction accuracy predominantly depends on the accuracy of the overall modelled mass and nacelle inertia, and provide equivalent tower dominated frequency predictions to the more complex flexible rotor model.
4. Model Updating
Foundation properties, , are estimated through a model updating approach, which optimises the compliance parameters of the rigid-rotor FE numerical model to achieve the best prediction of the modal behaviour identified from the field monitored data. The objective function shown in equation (29) is formed based on the sum of the normalised differences between the natural frequencies of the model fFE and those identified from the modal analysis of the structure fData. Mode shape comparisons, such as through modal assurance criteria (MAC), are not included in the objective function, due to the distortion of the closely spaced 1st tower mode shapes [35] and potential for tower-rotor coupling that is difficult to accurately simulate.
The OWT finite element model features a rigid-rotor-nacelle assembly, in which the rotor inertia about the axis of the turbine is neglected, representing idling conditions. The first three side-side and the 1stfore-aft tower mode are the four modes used to form the objective function and have been collectively identified in 64 of the available 87 idling datasets.
The foundation is modelled as an equivalent macroelement at seabed level, and the compliance terms are optimised during the model updating process. During this study, it was observed that model updating of the compliance components (C11, C12, and C22) provides faster and more robust optimisation towards global minima of the objective function, when compared to optimisation of the macroelement components directly (k11, k12, and k22). This observation is supported by Anderson [10], who concluded that the compliance components are better conditioned for optimisation than the stiffness components. It was therefore decided not to adopt any additional regularisation term in the objective function. Through summation of the objective function over the idling datasets in which the four target modes have been successfully identified, a single set of ensemble compliance components, , is optimised.
The MATLAB global optimisation routine, Pattern Search, is adopted to minimise the objective function. The allowable bounds of the compliance components are listed in Table 7. The bounds for the C11 and C22 terms are based on the typical values observed in the design models (shown in Table 5) but span approximately an order of magnitude above and below. The bounds for the coupling term C12 span a single order of magnitude to avoid solutions with unrealistically low coupling.
Three sets of initial estimates for the compliance parameters were trialled, which are the foundation parameters predicted by the API/DNVGL and PISA foundation models and a “high stiffness” set of parameters, as listed in Table 5. The “high stiffness” parameter set has been chosen to represent a stiff parameter set relative to both the PISA and API/DNVGL design values.
The average of the ensemble compliance components from model updating using the three initial parameter estimates, , is presented in Table 5. The three initial parameter estimates resulted in optimised parameters that differed by no more than 1.4%, suggesting a set of nearby minima identified from different initial estimates. Table 8 presents a comparison of the accuracy in finite element model prediction of the identified behaviour for the 64 idling datasets, when using the foundation design models or parameters. avg is the average of the normalised modal frequencies, is the average absolute percentage error of the individual dataset predicted frequencies, is the variance of eῶ, and are the average of the individual dataset modal assurance criteria, and S2 is the modal assurance criteria [22].
Across the three design model and foundation parameters, the 3rd tower fore-aft frequency was poorly predicted, which is attributed to highly sensitive blade-tower interaction of this mode.
While generally underpredicting stiffness, the PISA design model predicts the modal frequencies well, with around 0.5% average error for the 1st pair of tower modes and average errors of no more than 1.1% for the 2nd pair of tower modes. The MAC values between side-side identified and the PISA-predicted modes were above 0.9 indicating good prediction of the mode shape for all side-side modes. The close spacing of the 1st tower dominated modes is expected to result in a significant transformation of the resulting mode shapes. However, the high S2MAC values indicate that the observed modes stem from the subspace spanned by the two predicted mode shapes. The larger error in the predicted 2nd FA mode shape is also anticipated to occur due to the strong coupling of the blade and tower dynamics in the fore-aft direction.
Both the API/DNVGL and the Kallehave foundation models significantly underpredicted the identified natural frequencies and MAC values were typically marginally lower than for the PISA model, although still over 0.9. While the Kallehave model represented an improvement over the basic API/DNVGL model, it provided a substantially softer response prediction than the PISA model.
Predictions from all of the foundation models and parameters resulted in low variance in the frequency error, , indicating good robustness of the system model and modal identification data.
The parameters provided the lowest overall error in the predicted natural frequencies. However, mode shape accuracy was very similar to the PISA model. An example comparison of the observed and foundation finite element predicted mode shapes for a single dataset is presented in Figure 14, where good fidelity is observed at the sensor locations for all modes.
Figure 15 compares the effective kH and kM foundation stiffness at ground level against equivalent load eccentricity h, when using either design model foundation parameters or the parameters. The stiffness profiles reflect the predicted natural frequencies in Table 8, with the API/DNV design model grossly underpredicting kH and kM at all relevant load eccentricities. For values of h that are relevant to wind turbine dynamics, the response is predominantly dependent on kM, which is most heavily influenced by the C22 compliance term. This is demonstrated by the PISA effective stiffnesses, which overpredicts kM at low h and underpredicts at high h, resulting in a corresponding under- and overprediction of the 1st and 3rd tower natural frequencies, respectively, regardless of the overprediction of kH for all h. Similarly, whilst the API-Kallehave method indicates good prediction of kH, the underprediction of kM across all h results in underpredictions of the system natural frequencies. Figure 15 illustrates how good quality modal data for the 3rd and higher tower frequencies provides the information necessary to accurately identify the C11 parameter (the kH profile intercepts with the x-axis).
The overprediction of kH for the PISA model may have occurred due to the use of geotechnical site investigation data tailored towards an API/DNVGL design. Higher quality measurements of small strain shear stiffness (G0) variations in the field would likely improve the accuracy of the model prediction.
5. Conclusion
The paper demonstrates an approach to estimate foundation stiffness, based on finite element model updating against field monitored modal behaviour. Output only modal analysis, using the N4SID algorithm, provides stable modal estimates across the range of frequencies spanned by the first three tower dominated modes. This was repeated for 87 idling datasets over a four day period. The identified frequencies were shown to vary smoothly over time in accordance with the expected change of the environmental conditions. Complex normalisation of the mode shapes indicates a reduction in distortion effects associated with nonproportional damping.
Geotechnical site investigation data are used to predict the seabed level equivalent macroelement stiffness components for the API/DNVGL approach, the API/DNVGL approach with a Kallehave modification for sand initial modulus, and the PISA design framework. An approach is outlined for presenting the foundation stiffness, in terms of effective lateral and rotational stiffness against equivalent load eccentricity, which allows the relative stiffness at load eccentricities corresponding to relevant model behaviour to be compared.
Numerical Timoshenko beam finite element modelling of the structure, using an idling flexible rotor model based on the scaled NREL 5 MW baseline turbine, indicates that the decoupling of the blade and tower dynamics in the side-side direction allows the side-side tower dominated modal frequencies to be well predicted using a rigid-rotor-nacelle-assembly model. A model updating approach is applied to estimate the compliance components of the seabed level equivalent foundation macroelement, which indicates consistent convergence and independence of the initial compliance parameter estimates. Comparison of the model updated foundation parameters against the three design models indicates that the PISA model provides the most accurate prediction of the foundation stiffness and resulting turbine modal properties.
6. Limitations and Future Work
The presented model updating optimised foundation stiffness parameters and assumed characteristic design parameters for the monopile, tower, and nacelle. Future work should apply a suitable framework, such as through Bayesian inference methods, to investigate the influence of uncertainty in the relevant superstructure parameters and water added mass coefficient on the identified foundation parameters.
The study has investigated only idling datasets and where necessary has assumed blade mass and stiffness properties that are scaled from a benchmark turbine model. Future work should investigate the ability to identify foundation properties during more common operational periods, the accuracy of which will be critically dependent on accurate modelling of the mass and stiffness properties of the rotor-nacelle assembly.
It is evident that the foundation effective lateral and moment stiffnesses are complicatedly coupled to the compliance terms and the modal behaviour. Whilst Anderson [10] investigated the sensitivity of optimisation error to the individual compliance components, further work should expand on and demonstrate the sensitivity of OWT natural frequencies to the foundation effective lateral and moment stiffnesses, to demonstrate how the accuracy of the different components will influence loads and design.
Data Availability
The work presented in this article is based on commercially confidential data provided by Parkwind NV and is not available to be shared.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work presented in this paper was supported by the EPSRC Supergen Wind Hub programme (Grant no. EP/L014106/1) through their flexible funding scheme. The authors gratefully acknowledge the support of Parkwind NV for granting access to the monitoring data and OWI Lab for providing data files and design documents. Byrne is supported by the Royal Academy of Engineering under the Research Chairs and Senior Research Fellowships scheme.