Comparison of measured and simulated fatigue loads on a multi-megawatt wind turbine

The Mann turbulence model is widely used in the design and certification of multi-megawatt wind turbines. However, these turbines operate in a region of the atmosphere where the model's assumptions are violated. One of the most significant assumptions is that of neutral stability conditions, which raises concerns about the model's accuracy for load simulations. To investigate this, we compare fatigue loads measured on a 15 MW wind turbine to simulations performed using an aeroelastic solver. The inflow was characterized using data from a meteorological mast equipped with sonic and cup anemometers. The turbulence model was fitted to measurements of auto-spectra under varying wind speeds and stability conditions, while the vertical profile of wind speed was represented by a power law. The resulting wind fields were then used as input to the aeroelastic simulations.

We first present a comparison of measured fatigue loads on the tower and blades across different atmospheric stability regimes. On average, the fore-aft loads at the bottom of the tower were 98 % higher under unstable atmospheric conditions as compared to stable conditions. On the other hand, the difference in flapwise blade loads between unstable and stable conditions was 20 % below rated wind speed and −2 % at higher wind speeds. These results underscore the importance of accounting for atmospheric stability in wind turbine siting and load verification campaigns. A subsequent comparison of measurements and simulations revealed that measured and simulated loads tend to fall within 3 standard deviations of each other, even under non-neutral conditions. However, the simulated fatigue loads on the tower were overestimated by a margin of 5 standard deviations under some stable conditions, likely due to incorrect predictions of the spectral coherence made by the turbulence model. Shear extrapolation based on the power law might also lead to overestimation of blade loads in the simulations.

These results indicate that, despite its simplifying assumptions, the Mann model, when fitted to measurements of turbulence auto-spectra, does not introduce significant errors in fatigue load simulations for solitary multi-megawatt wind turbines.

1 Introduction

The Wind Turbines – Part 1: Design Requirements of the International Electrotechnical Commission (IEC) (IEC, 2019) recommends the Mann (Mann, 1994) turbulence model and IEC Kaimal (Kaimal et al., 1972) spectra with exponential decay coherence (Davenport, 1961) as inputs for aeroelastic simulations of wind turbines. These models are designed to capture turbulence in the surface layer of the atmosphere. However, modern wind turbines with rotor diameters greater than 200 m and hub heights exceeding 100 m often operate outside the surface layer. The assumptions behind the aforementioned models are no longer valid in this region of the atmosphere (Veers et al., 2023, Sect. 3.2, p. 1083). For instance, the models assume neutral stability conditions, homogeneous turbulence in all directions of space and Gaussian statistics at the smallest scales. On the contrary, atmospheric turbulence is often thermally generated (Obukhov, 1962), is inhomogeneous in the vertical direction and follows non-Gaussian statistics at the smallest scales (Wyngaard, 2010, p. 156). Thus, it is uncertain whether the Mann and Kaimal models can provide accurate estimates of the dynamic loads on multi-megawatt wind turbines.

The aim of this study is to gauge the accuracy of load simulations where the input turbulence field is described by the Mann model, by comparing with measurements of fatigue loads from a 15 MW prototype wind turbine.

The impact of inflow turbulence on the fatigue loads experienced by a wind turbine has been well documented in the literature. Robertson et al. (2019) performed a sensitivity study to evaluate the influence of various inflow parameters such as turbulence intensity, length scales, vertical shear and veer. They found that the turbulence and length scale of the streamwise wind component, along with vertical shear, played a primary role in determining the fatigue loads on a turbine. Sathe et al. (2013) theoretically studied the impact of atmospheric stability on fatigue loads. They showed that the tower and rotor loads were significantly higher during unstable conditions compared to the loads under stable stratification. Nybø et al. (2021a) simulated the fatigue loads on a 15 MW wind turbine using different inflow specifications such as turbulence boxes generated using the Mann and Kaimal models and wind fields from large-eddy simulations (LESs). The simulated loads showed considerable differences even when the three methods had the same one-point statistics. The variability was attributed to the difference in the spatial structure of turbulence, but given the lack of loads measurements, it was not possible to determine which inflow models were the most accurate.

Concurrently, the Mann model has been modified to capture phenomena that were lacking in the original formulation. These include atmospheric stability (Chougule et al., 2018), dependence on time (de Maré and Mann, 2016), mesoscale turbulence (Syed and Mann, 2024) and inhomogeneity (Guo et al., 2026), although the model of Chougule et al. (2018) is not numerically stable for convective conditions. On the other hand, non-Gaussian turbulence and effects of the Coriolis force have not been included in the model as of yet. However, given the absence of measurements of loads from an operational wind turbine to serve as the ground truth, it is unclear whether these phenomena lead to significant errors in aeroelastic simulations. Resolving this uncertainty requires a so-called “one-to-one” comparison, where measured loads are compared to simulated loads. The actual inflow experienced by the turbine is replicated in the simulation by using measurements of turbulence from a met-mast or a lidar. In this regard, Brown et al. (2024) performed a one-to-one comparison on a 2.8 MW turbine. Inflow measurements from an array of sonic anemometers, wind vanes and a nacelle lidar were assimilated into the simulation environment using a few different methods (Jonkman and Buhl, 2006; Rinker, 2022). This study showcased the value of one-to-one comparisons while also highlighting the associated difficulties. The flapwise blade and fore-aft tower loads from the simulations showed a significant bias compared to the measurements, due to inaccuracies in the aerodynamic modelling and the absence of certain controller features in the simulation environment. Indeed, Doubrawa et al. (2024) showed that the bias can be reduced to a large extent by incorporating some features of an industrial controller such as rotor position control. Although in most studies, it is difficult to replicate the turbine in a simulation as the blade design and controller architecture are confidential. However, the work of Dimitrov and Natarajan (2017) did not suffer from these drawbacks. They compared inflow assimilation using met-masts to lidar-based methods and found that the uncertainties in the load estimates from lidar-based and met-mast-based flow assimilation were similar and occasionally lower when using the lidar. Moreover, the simulated loads agreed with the measurements at different wind speeds and turbulence intensities. However, this study considered an older turbine variant with all measurements taken below 100 m in the atmosphere. Thus, the validity of the Mann model for load simulations of multi-megawatt wind turbines remained an open question.

We present fatigue loads measured on a 15 MW prototype under varying atmospheric conditions and compare them to aeroelastic simulations with the same controller architecture, and aerodynamic and structural properties as the real turbine. The inflow is assimilated by fitting the Mann model to measurements of turbulence auto-spectra from a 155 m tall met-mast. To the authors' knowledge, the relationship between fatigue loads on a wind turbine operating in the field and atmospheric stability has previously not been documented in the literature. Moreover, we present full-scale measurements in the natural atmosphere which form the basis for load validation. Full-scale measurements are important in this context as they capture features of the atmosphere that are not modelled by simulations.

The paper is organized as follows: Sect. 2 presents the measurement campaign, the method of inflow assimilation and a brief overview of the simulation environment. The results are described in Sect. 3, starting with the dependence of fatigue loads on atmospheric stability, followed by one-to-one comparisons of blade and tower loads. Finally, the conclusions of this study are reflected upon in Sect. 4.

2 Methods

2.1 Measurement campaign

The purpose of the measurement campaign was, among others, power curve verification and load validation of the Vestas V236 prototype wind turbine. The turbine has a rotor diameter of 236 m and a hub height of 163 m, with a rated power output of 15 MW. It is an offshore wind turbine, and several variants are currently under installation in wind farms in the North Sea and the Baltic Sea.

The campaign lasted for 2 years, from 2023 to 2025, and took place at the Test Center Østerild. It is located in north Jutland, Denmark, about 50 km from the coast. We refer the reader to Peña (2019) for a detailed description of the site. The prototype turbine was equipped with many sensors that measured bending moments on various components, the power output of the turbine, pitch angle of all the blades and rotation speed, among other operational data. An overview of the data signals used in this study is shown in Table 1. An IEC-complaint meteorological mast was also present, about 660 m to the west of the prototype. Since the dominant wind direction was westerly, the mast was often present upstream of the turbine. It was equipped with a sonic anemometer at 155 m; cup anemometers at 45 m, 104 and 163 m; and wind vanes at 45 and 159 m above ground level.

Table 1 An overview of the signals analysed in this study and the devices that measure them.

2.2 Processing and binning of data

The steps involved in the curation of the data are described below, with the data availability after each step shown in Table 2. Note that the goal of data curation was to obtain turbulence measurements appropriate for a one-to-one comparison.

Table 2 The data processing steps and the resulting availability. Note that the data availability in each subsequent row is expressed as a percentage of the total data available between the start and end of the campaign.

The data were sampled at a frequency of 5 Hz and partitioned into periods of 10 min. This sampling frequency is chosen as turbulent fluctuations that are faster than the Nyquist frequency of 2.5 Hz are not deemed to be important from a fatigue loads perspective. Only those periods were selected for further analysis, where the 10 min mean wind direction $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ measured by the wind vane at 163 m was between 225° and 315°. Note that $\stackrel{\mathrm{‾}}{\cdot }$ refers to time averaging, and the wind directions are expressed in meteorological convention, with 0° being wind coming from the north. The westerly sector is selected because winds from this direction can be assumed to be homogenous and free from the wake of neighbouring turbines (Peña, 2019). Subsequently, the 10 min mean wind speed ($\stackrel{\mathrm{‾}}{u}$) and direction measured by the sonic anemometer at 155 m were compared to the cup and vane at 163 and 159 m, respectively. The purpose of this comparison was to identify inaccuracies in the instruments and remove the corresponding periods from the dataset. As seen in Fig. 1, the instruments agreed in the measurements of the mean wind speed and direction. Although biases of 0.3 m s⁻¹ and 1° are detected in (a) and (b), respectively, this can be attributed to the sonic anemometer being placed 8 and 4 m lower than the cup and vane, respectively. Next, the 10 min periods where the mean wind speed was outside the turbine operating region were removed. This was followed by calculating the Monin–Obukhov length (L_mo) for each of the remaining 10 min periods:

$$\begin{array}{}\text{(1)}& L_{\text{mo}}={\displaystyle \frac{-u_{*}^{\mathrm{3}}\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v}}}}{\mathit{\kappa }g\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}}},\end{array}$$

where θ_v is the virtual potential temperature, κ is the von Kármán constant, w^′ denotes the fluctuating vertical wind-velocity component and $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ is the vertical flux of virtual potential temperature. Moreover, the friction velocity, u_* is given by

$$\begin{array}{}\text{(2)}& u_{*}=\sqrt{-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}.\end{array}$$

Note that we used the sonic temperature (T_s) in Eq. (1), as $\stackrel{\mathrm{‾}}{T_{\mathrm{s}}}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{T}_{\mathrm{s}}^{\prime }}$ are considered to be acceptable approximations to $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v}}}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$, respectively (Schotanus et al., 1983).

Figure 1 (a) The 10 min mean wind speed measured by the sonic anemometer at a height of 155 m compared to measurements made by the cup anemometer at a height of 163 m, when the mean wind direction was between 225° and 315°. (b) The 10 min mean wind direction (in meteorological coordinates) measured by the sonic anemometer compared to measurements made by the wind vane at a height of 159 m. Note that the red line has a slope equal to 1.

The 10 min data were subsequently binned according to the mean wind speed measured by the sonic anemometer at 155 m and the Monin–Obukhov length based on the stability classifications suggested by Gryning et al. (2007) and shown in Table 3. The size of the mean wind speed bins was 2 m s⁻¹, with the bin centres being odd numbers from 5–21. Table 3 is only applicable within the surface layer where the friction velocity and temperature flux do not change with height. While this assumption may not be strictly valid for the measurement height of 155 m, the reader is cautioned against using Table 3 outside the surface layer.

Table 3 Classifications of atmospheric stability based on the Monin–Obukhov length.

The probability distribution of L_mo conditioned on $\stackrel{\mathrm{‾}}{u}$ is shown in Fig. 2. We found stable conditions to be most frequent at low wind speeds, contrary to the observations of Peña (2019), wherein unstable conditions were dominant at wind speeds below 11 m s⁻¹. This observation is likely a sampling bias. The instruments briefly stopped recording data in the summer of 2024 (as noted in row 2 of Table 2) when unstable conditions are more frequent. Neutral and near-neutral conditions were dominant at mean wind speeds above 14 m s⁻¹, which is likely a result of mechanically generated turbulence being relatively stronger than the effects of buoyancy.

Figure 2 The distribution of atmospheric stability quantified by the Monin–Obukhov length (L_mo) conditioned on the 10 min mean wind speed ($\stackrel{\mathrm{‾}}{u}$).

In the subsequent processing step, we aimed to find the longest period of contiguous measurements within each bin. Table 4 shows the capture matrix indicating the number of contiguous 10 min periods (N_bin) found for each bin. The turbulence measurements (namely the auto-spectra) were averaged over N_bin periods and fitted to the Mann model. Thus, we used Bartlett's method (Bartlett, 1948) to compute the auto-spectra. If N_bin<4, the uncertainties in the measurements were deemed to be too high and a fit with the Mann model could not be obtained. As a consequence, a one-to-one comparison could not be carried out in the corresponding bins. In other words, each bin required at least 40 min of continuous data. This condition was satisfied in 39 bins.

The choice to limit the analysis to periods with contiguous data was motivated by the observation that wind speed and stability did not fully describe turbulence at these heights. Indeed, when computing the spectra for all 10 min periods within a given bin, we found the spread to be large, encompassing 3 orders of magnitude, from 10⁻³ to 1 m²s⁻². Additionally, as described later, spectral fitting was the preferred method for inflow assimilation. Thus, by attempting to fit the Mann model to the bin-averaged spectra, we risked lumping together distinct atmospheric conditions which would add significant uncertainty to the comparisons of loads. On the other hand, by selecting a contiguous period of measurements with minimal spread in the individual spectra, the resultant uncertainty was lower and easier to quantify as it was purely statistical in nature.

2.3 Assimilation of inflow

The goal of inflow assimilation is to achieve a faithful representation of the measured inflow in the simulation environment. To replicate the turbulence in the simulations, we use spectral fitting. Herein, the Mann model is fitted to the auto-spectra of the along-wind u, cross-wind v and vertical w wind components as well as the real part of the uw cross-spectrum to obtain the three model parameters: L, $\mathit{\alpha }{\mathit{\epsilon }}^{\mathrm{2}/\mathrm{3}}$ and Γ for each stability and wind speed bin. These parameters are then used to generate a turbulence box which is provided as input to the aeroelastic solver (Mann, 1998).

Previous one-to-one comparisons have used different inflow assimilation methods such as scaling a turbulence box with the measured variances in the three components or by constraining a turbulence box on a set of measurements. However, the former has been shown to possess a high level of uncertainty (Liew and Larsen, 2022), while the latter is more suitable for measurements from a nacelle lidar (Rinker, 2022). Since measurements of the auto-spectra contain the length scale and variance, both of which have a significant influence on fatigue loads, fitting the turbulence model to the spectra can achieve accurate assimilation of the inflow for the purposes of load validation (Nybø et al., 2021b).

Table 4 The number of contiguous 10 min periods (N_bin) per mean wind speed and stability bin used to compute the auto-spectra of turbulence.

The fitting procedure is as follows. Measurements of auto-spectra are computed by

$$\begin{array}{}\text{(3)}& S_i\left(f\right)={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{N_{\text{bin}}T}}\sum _{k=\mathrm{1}}^{N_{\text{bin}}}{\left|{\mathrm{u}}_i^k(f,T)\right|}^{\mathrm{2}}\end{array}$$

where

$$\begin{array}{}\text{(4)}& u_i^k(f,T)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{T}{\int }}{u}_i^k\left(t\right)\exp (-\mathrm{2}\mathit{\pi }ift)\mathrm{d}t,\end{array}$$

is the Fourier transform of the ith component such that $(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}})=(u,v,w)$ and T=600 s. Note that measurements of the three wind components at a sampling frequency of 5 Hz are obtained from the sonic anemometer. The spectra are further block-averaged over logarithmic bins N_l to reduce the level of noise at the high frequency end of the spectrum. The block-averaging is done such that there are 12 points per decade, which means that N_l increases with frequency. Hence, the statistical uncertainty in the measurements σ(〈S_i〉) is quantified by

$$\begin{array}{}\text{(5)}& \mathit{\sigma }(\langle S_i\rangle )={\displaystyle \frac{\langle S_i\rangle }{\sqrt{N_{\text{bin}}{N}_{\mathrm{l}}}}}.\end{array}$$

Then, the model parameters that minimize the error ξ² are obtained by an optimization algorithm. ξ² is defined as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\mathit{\xi }}^{\mathrm{2}}\left(L,\mathrm{\Gamma },\mathit{\alpha }{\mathit{\epsilon }}^{\frac{\mathrm{2}}{\mathrm{3}}}\right)={\sum }_{i=\mathrm{1}}^{\mathrm{3}}{\sum }_{n=\mathrm{1}}^N{\left[\log \left(k_{\mathrm{1}}^n{F}_i\right)-\log \left(k_{\mathrm{1}}^n{F}_{i,m}\right)\right]}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{n=\mathrm{1}}^N{\left[\log \left(k_{\mathrm{1}}^n\mathrm{\Re }\left({\mathit{\chi }}_{uw}\right)\right)-\log \left(k_{\mathrm{1}}^n\mathrm{\Re }\left({\mathit{\chi }}_{uw,m}\right)\right)\right]}^{\mathrm{2}},\end{array}\end{array}$$

where the subscript m refers to the estimated spectra and χ_uw is the uw cross-spectrum. On the other hand, the modelled spectra are given by

$$\begin{array}{}\text{(7)}& F_i\left(k_{\mathrm{1}};L,\mathrm{\Gamma },\mathit{\alpha }{\mathit{\epsilon }}^{\frac{\mathrm{2}}{\mathrm{3}}}\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}{\mathrm{\Phi }}_{ii}\left(\mathit{k};L,\mathrm{\Gamma },\mathit{\alpha }{\mathit{\epsilon }}^{\frac{\mathrm{2}}{\mathrm{3}}}\right)\mathrm{d}{k}_{\mathrm{2}}\mathrm{d}{k}_{\mathrm{3}}.\end{array}$$

Here, Φ(k) is the spectral tensor derived by Mann (1994) defined over wave number space: $\mathit{k}=(k_{\mathrm{1}},k_{\mathrm{2}},k_{\mathrm{3}})$. Note that Taylor's hypothesis is used to move between frequency and wave number space. Moreover, the coherence, ${\mathit{\gamma }}_{uu}^{\mathrm{2}}$, is defined as

$$\begin{array}{}\text{(8)}& {\mathit{\gamma }}_{uu}^{\mathrm{2}}(f,{\mathrm{\Delta }}_y,{\mathrm{\Delta }}_z)={\displaystyle \frac{\left|{\mathit{\chi }}_{uu}\right(f,{\mathrm{\Delta }}_y,{\mathrm{\Delta }}_z){|}^{\mathrm{2}}}{S_u(f,\mathrm{0},\mathrm{0})S_u(f,{\mathrm{\Delta }}_y,{\mathrm{\Delta }}_z)}},\end{array}$$

where Δ_y and Δ_z are the lateral and vertical separations, respectively, and χ_uu is the cross-spectrum between $u(t,\mathrm{0},\mathrm{0})$ and $u(t,{\mathrm{\Delta }}_y,{\mathrm{\Delta }}_z)$.

Figure 3 shows the estimated spectra under different stability conditions along with the fit derived from Eq. (6) at a mean wind speed of 9 m s⁻¹. Note that the area under the spectra of a given component is equal to the variance of that component, while the location of the spectral peak is related to the length scale. As expected, unstable conditions show much more variance (or turbulence intensity) and larger length scales as compared to stable conditions, whereas neutral conditions are intermediate. The Mann model fits the spectra well in the given frequency range.

Figure 3 The estimated pre-multiplied auto-spectra of the u, v and w components and the real part of the uw cross-spectrum under varying stability conditions, along with a mean wind speed of 9 m s⁻¹.

The model parameters obtained in each stability and wind speed bin are shown in Fig. 4. The parameter $\mathit{\alpha }{\mathit{\epsilon }}^{\mathrm{2}/\mathrm{3}}$ is proportional to the square of the mean wind speed (Syed and Mann, 2024) and hence $\sqrt{\mathit{\alpha }{\mathit{\epsilon }}^{\mathrm{2}/\mathrm{3}}}$ is observed to vary linearly with $\stackrel{\mathrm{‾}}{u}$. Furthermore, $\mathit{\alpha }{\mathit{\epsilon }}^{\mathrm{2}/\mathrm{3}}$ is also related to the variance in the wind speed (Kelly, 2018). Thus, Fig. 4a shows that higher wind speeds are associated with more turbulence. Length scales as large as 180 m are observed in certain unstable conditions, while stable conditions usually show length scales lower than 50 m. The anisotropy parameter Γ appears to be lower for unstable conditions as compared to stable conditions at low wind speeds. Γ is related to the ratio of the variances in the three wind components, and Γ=0 when ${\mathit{\sigma }}_u={\mathit{\sigma }}_v={\mathit{\sigma }}_w$. Thus, convective conditions show lower values of Γ, as the fluctuations in v and w tend be to stronger than under stable conditions. The variation in Γ decreases with wind speed, probably due to the decrease in the data availability of non-neutral stabilities at higher wind speeds. Note that Γ might also depend on surface roughness and boundary layer depth, which are outside the scope of this study but can be investigated in the future.

Figure 4 Parameters of the Mann turbulence model as functions of wind speed and stability obtained by fitting the model to the measurements shown in Table 4.

Given the importance of vertical shear and veer on fatigue loads (Robertson et al., 2019), we also attempted to replicate the measured variation of mean wind speed and direction with height in the simulation setup. Data from three cup anemometers present at heights 45, 105 and 163 m, and two wind vanes present at 45 and 159 m were used to obtain the mean vertical profiles of the along-wind component and the wind direction $\stackrel{\mathrm{‾}}{u}\left(z\right)$ and $\stackrel{\mathrm{‾}}{\mathit{\theta }}\left(z\right)$. Next, the shear parameter α_s was extracted by fitting the power law profile to the measurements:

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{‾}}{u}\left(z\right)=\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{h}}\right){\left({\displaystyle \frac{z}{z_{\mathrm{h}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{s}}},\end{array}$$

where z_h is the hub height. Here, the shear above the rotor is assumed to follow the same power profile as measured below the rotor. Note that Eq. (9) provided a good fit to the measurements as the coefficient of determination (R²) was always greater than 0.92 in all 39 bins. The veer profile, on the other hand, is assumed to be a linear function of height:

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{\mathit{\theta }}\left(z\right)=\mathit{\beta }(z-z_{\mathrm{h}}).\end{array}$$

β is obtained by fitting the measured profile of $\stackrel{\mathrm{‾}}{\mathit{\theta }}\left(z\right)$ to Eq. (10). The values of α_s and β for all stability and wind speed bins are shown in Fig. 5. A larger variation with height in the wind speed and direction is observed during stable conditions as compared to unstable conditions, as more turbulent mixing occurs in the latter, which leads to more uniformity in the mean wind speed and direction.

Figure 5 The shear and veer parameters as functions of wind speed and stability obtained by fitting Eqs. (9) and (10) to measurements of the vertical profiles of mean wind speed and direction.

Thus, the wind box u_p provided to the aeroelastic solver is

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle }{\mathit{u}}_{\mathit{p}}(x,y,z)=& {\displaystyle }{\mathit{u}}_{\mathbf{m}}(x,y,z;\mathit{\alpha }{\mathit{\epsilon }}^{\mathrm{2}/\mathrm{3}},L,\mathrm{\Gamma })+{\mathit{e}}_x\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{h}}\right){\left({\displaystyle \frac{z}{z_{\mathrm{h}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{s}}}\\ {\displaystyle }& {\displaystyle }+{\mathit{e}}_y\stackrel{\mathrm{‾}}{u}\left(z\right)\tan \left(\mathit{\beta }\right(z-z_{\mathrm{h}}\left)\right).\end{array}\end{array}$$

Here, u_m is the wind box of turbulent fluctuations with zero mean generated using the Hipersim library (Dimitrov et al., 2024) and e_x, e_y are the unit vectors in the along-wind and cross-wind directions.

The turbulence boxes had a temporal resolution of Δt=0.15 s, hence the spatial resolution in the along-wind direction Δx was $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{h}}\right)\mathrm{\Delta }t$. Given a total simulation time of 600 s, the number of grid points in the x direction were 8094. On the other hand, the vertical and lateral grid resolution was 3.8 m, with 64 grid points in either direction. These values were recommended by Liew and Larsen (2022) for the convergence of fatigue load estimates. They also suggest using at least seven turbulence seeds to reach a standard error of less than 5 % in the damage equivalent moment on the blades and the tower. Thus, in this study, a total of $\mathrm{7}\times \mathrm{39}=\mathrm{273}$ simulations were performed.

2.4 A note on the aeroelastic solver

The dynamic response of the turbine to a given wind input is simulated in the Vestas turbine simulator (VTS). It is a state-of-the-art blade element and momentum theory (BEM) solver derived from Flex5 (Øye, 1996) and based on the mode shape approach using generalized coordinates, masses, stiffnesses and forces. The solver is non-linear as mass, damping and stiffness matrices are recalculated at every time step. It includes tip loss correction and a dynamic stall model, and is linked to the controller via a dynamic-link library (DLL). Since the system of equations is not stiff, a very efficient equation solver may be used (Runge–Kutta–Nyström with Cholesky decomposition), with a relatively large constant time step.

2.5 Calculation of fatigue damage and quantification of load differences

Fatigue loads are often characterized by the 1 Hz damage equivalent moment D_M, which is based on the Palmgren–Miner damage rule (Burton et al., 2011). Consider S(t) to be the time-dependent bending moment on a particular component. S(t) is first decomposed into a set of single amplitude load ranges S_i, and the number of cycles for each load range n_i are calculated using rain-flow counting (Rychlik, 1987). Moreover, the number of cycles to failure N_i at a certain load range are found using the Wöhler S–N curve:

$$\begin{array}{}\text{(12)}& N_i={\left({\displaystyle \frac{S_{\mathrm{0}}}{S_i}}\right)}^m,\end{array}$$

where m is the Wöhler exponent of the material and S₀ is the y intercept of the S–N curve. Note that in this study, m is taken as 4 for steel and 10 for glass fibre composite (Slot et al., 2019). The fatigue damage suffered by the material at particular load range is given by

$$\begin{array}{}\text{(13)}& D_i={\displaystyle \frac{n_i}{N_i}}.\end{array}$$

Thus, the total damage according to Palmgren–Miner's rule is

$$\begin{array}{}\text{(14)}& D=\sum D_i=\sum {\displaystyle \frac{n_i}{N_i}}.\end{array}$$

Finally, D_M is defined as the single amplitude time-varying load applied on the component for n_eq cycles such that it results in the same level of accumulated damage as S(t):

$$\begin{array}{}\text{(15)}& D_{\mathrm{M}}\left(S\right)={\left({\displaystyle \frac{\sum n_i{S}_i^m}{n_{\text{eq}}}}\right)}^{\frac{\mathrm{1}}{m}},\end{array}$$

where n_eq is 600, which implies that D_M(S) is applied at a frequency of 1 Hz for a duration of 10 min.

We inform the reader that all the results presented in the next section have been anonymized. This was achieved by scaling the absolute value of a given quantity (say X) as follows:

$$\begin{array}{}\text{(16)}& \stackrel{\mathrm{̃}}{X}={\displaystyle \frac{X}{X_{\text{max}}-X_{\text{min}}}},\end{array}$$

where X_min and X_max are the minimum and maximum values of X, respectively. Note that the minimum and maximum are the same for a given load channel, and are computed over all stabilities and mean wind speeds. Thus, all loads from a particular channel, irrespective of whether they are obtained from measurements or simulations, are scaled with the same values of X_min and X_max. As a result, Eq. (16) is a rescaling operation that preserves all qualitative differences between the data points.

The difference between damage equivalent moments in the simulations and measurements are quantified using the z score Z defined as follows:

$$\begin{array}{}\text{(17)}& Z={\displaystyle \frac{\stackrel{\mathrm{‾}}{A}-\stackrel{\mathrm{‾}}{B}}{\sqrt{{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{A}}^{\mathrm{2}}+{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{B}}^{\mathrm{2}}}}},\end{array}$$

where

$$\begin{array}{}\text{(18)}& {\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{A}}={\displaystyle \frac{{\mathit{\sigma }}_A}{N_A}}.\end{array}$$

A and B refer to the measured and simulated quantities respectively. N_A is the number of 10 min samples whose values can be found in Table 4, while N_B=7, which corresponds to the number of seeds used in each wind speed and stability bin. Note that when Z=1, $\stackrel{\mathrm{‾}}{A}$ and $\stackrel{\mathrm{‾}}{B}$ are 1σ apart, where

$$\begin{array}{}\text{(19)}& \mathit{\sigma }=\sqrt{{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{A}}^{\mathrm{2}}+{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{B}}^{\mathrm{2}}}.\end{array}$$

Furthermore, a positive value indicates that fatigue loads are overestimated in the simulation, while a negative value indicates an underestimation.

3 Results

3.1 Effect of atmospheric stability on fatigue loads

The mean damage equivalent moments are computed from the measurements of bending moments at the tower bottom and blade root, and are shown in Figs. 6 and 8 as functions of wind speed and atmospheric stability. Note that the fore-aft direction is parallel to the mean wind direction and side-side is perpendicular to it, when the yaw-misalignment is zero. Moreover, Fig. 8 shows the damage equivalent loads on one of the blades, although similar results were obtained for all three. Lastly, the error bar indicates the standard deviation of the damage equivalent moment, which is also scaled according to Eq. (16).

Figure 6 The mean damage equivalent moments on the tower bottom in the (a) fore-aft and (b) side-side directions as functions of wind speed and stability. Note that the data have been anonymized using Eq. (16), and the error bars indicate the standard deviation. Moreover, (a, b) are scaled with different values of X_min and X_max.

Figure 6a shows that the damage equivalent moments at the tower bottom in the fore-aft direction are significantly influenced by atmospheric stability. For instance, in the 11 m s⁻¹ wind speed bin, which is around the turbine's rated wind speed, the fatigue loads are twice as high under unstable and neutral conditions as compared to stable conditions. Averaging over the wind speed bins where data are simultaneously available for unstable, neutral and stable conditions, the damage equivalent moments are 98 % and 65 % higher under unstable and neutral conditions as compared to stable ones. Moreover, the standard deviation of the damage equivalent moment is also 1.5 times higher under unstable and neutral conditions. On the other hand, the fatigue loads in the side-side direction are not as strongly influenced by stability. The difference in damage equivalent loads in this direction between unstable and stable conditions is 17 %. Similar trends are observed in the damage equivalent moments at the tower top.

As seen in Figs. 3, 5 and 7, unstable conditions are generally associated with higher turbulence intensities (TI_u) and larger length scales but lower shear in the vertical wind profile. Thus, the results presented herewith show that tower loads are mainly driven by turbulence and not shear. Indeed, it is expected that a three-bladed rotor will average out the effects of load variation due to the wind profile. The TI_u, on the other hand, quantifies the fluctuations of a given wind component around its mean, and consequently, a larger TI_u of the along-wind component can lead to larger variations in tower loads. However, the TI_u alone is not a sufficient proxy for fatigue loads. The length scale describes the size of the eddies, which drive the turbulent fluctuations. Eddies that are comparable to the rotor diameter cause variations in the dynamic loads, while eddies that are smaller in size are filtered out. Indeed, these results are not surprising given the insights from buffeting theory for bridges and towers. Davenport (1964) and Scanlan (1978) showed that for a linear approximation, the buffeting response is proportional to the turbulence intensity. Hence, the TI_u and length scale together are responsible for fatigue loads on the tower being 98 % higher under unstable conditions as compared to stable conditions. We also note that the damage equivalent moments under neutral conditions do not necessarily represent the mean loads at any given wind speed. They appear to be closer to unstable conditions than to stable ones. This is because the TI and length scale of the load-driving u component under neutral conditions was closer to the values under unstable conditions.

Figure 7 The measured turbulence intensity TI_u of the along-wind component for each wind speed and stability bin at a height of 155 m. Note that TI_u is computed as the average over all 10 min periods in a given bin.

Figure 8 The mean damage equivalent moments at the root of blade 1 in the (a) flapwise and (b) edgewise directions as functions of wind speed at 155 m and stability.

The damage equivalent loads on the blade in the flapwise direction are impacted to a lower degree by atmospheric stability. At wind speeds below rated (the rated wind speed is approximately 11 m s⁻¹), unstable and neutral conditions cause approximately 20 % more fatigue damage as compared to stable conditions. However, at wind speeds above rated, the damage equivalent loads under neutral and stable conditions differ by −2 %. These observations are caused by the opposing effects of turbulence and shear, both of which affect the fatigue loads on a wind turbine blade. While unstable conditions are more turbulent, the mean wind speed does not vary appreciably with height. Stable conditions show converse behaviour such that the two drivers of fatigue loads (turbulence and shear) tend to nullify each other in Fig. 8. Moreover, as pitch control is activated at wind speeds above rated, the blade pitches out of the mean wind. Hence, the flapwise loads become less sensitive to turbulence and shear becomes relatively dominant. This causes stable conditions to be slightly more detrimental at high wind speeds. In the case of the edgewise loads on the blade, as they are mainly driven by gravity, no significant dependence on stability or turbulence is observed. Thus, we exclude them from the one-to-one comparisons presented in the next section. The side-side loads on the tower bottom are also excluded as they are not only dependent on the inflow but also on a range of different factors such as mass and pitch imbalances on the blades. As it is difficult to accurately quantify these imbalances, they cannot be replicated in a simulation. Hence, attributing any disagreement between measured and simulated loads to differences in the inflow cannot be achieved with a sufficient degree of confidence. On the other hand, the tower fore-aft and blade flapwise loads are significantly more sensitive to the inflow. Thus, our one-to-one comparison is limited to these two load channels.

3.2 One-to-one comparisons

Figure 9 shows the one-to-one comparison of damage equivalent moments on the bottom on the tower in the fore-aft direction, while Table 5 presents the z scores. They indicate that the tower-bottom damage equivalent moments are overestimated when inflow turbulence is modelled with the Mann model, and the model parameters are derived by fitting the model to measurements of the auto-spectra. In most instances, the overestimation is fairly low and the loads differ by less than 3σ. However, at certain wind speeds under near-neutral stable, stable, very stable and unstable conditions, the loads are overestimated by more than 5σ. We found that in the bins with near-neutral, stable and very stable conditions, the vertical coherence in the u component predicted by the model was higher than the measurements. Using data from the sonic anemometer at 155 m, and the cup anemometer and wind vane at 45 m, we found measurements of ${\mathit{\gamma }}_{uu}^{\mathrm{2}}(f,\mathrm{0},\mathrm{110})$, which were compared to the coherences from the turbulence model. In the frequencies between 10⁻³ and 10⁻² Hz, the model overestimated the coherence by up to 174 %. We suspect that the lateral coherence is also overestimated in these cases. Indeed, Nybø et al. (2021a) showed that coherence in the frequency range between 10⁻³ and 10⁻² Hz has a significant impact on the fatigue loads experienced by the towers of multi-megawatt turbines. Regarding the case with unstable conditions at 11 m s⁻¹, the measurements of the auto-spectra show that the variance in the w component is higher than that in the u and v components at frequencies below 10⁻² Hz. Such observations have been made before under strongly convective conditions (Syed, 2024) at an altitude of 200 m. However, the Mann model is based on the assumption that ${\mathit{\sigma }}_u>{\mathit{\sigma }}_v>{\mathit{\sigma }}_w$, which is true in the case of shear-driven turbulence but invalid for some unstable conditions, especially for altitudes above 100 m. Thus, fitting the model to the observations of the auto-spectra results in a relatively poor fit, which we believe causes the large mismatch between simulated and measured loads on the tower (and the blade). The one-to-one comparison of tower-top moments in the fore-aft direction were found to be identical to the tower bottom.

Figure 9 The damage equivalent moment at the bottom of the tower in the fore-aft direction, as measured on the prototype (in black) and simulated (in red) in the aeroelastic solver by assimilating the measurements of inflow via spectral fitting with the Mann turbulence model, along with a power law profile for the vertical variation in the wind speed.

Table 5 The z scores for the damage equivalent moments on the bottom of the tower in the fore-aft direction.

The one-to-one comparison of the damage equivalent moments on the blade in the flapwise direction is shown in Fig. 10, while the z scores are displayed in Table 6. In this case as well, the difference in loads is often less than 3σ. Although the measured loads in Fig. 10d appear to be higher than the simulated loads at 17 m s⁻¹ and higher, the absolute difference being less than 3σ means that the discrepancy can be better explained as a statistical error. However, at mean wind speeds of 9, 11 and 13 m s⁻¹, the damage equivalent moments under very stable and stable conditions are overestimated by more than 11σ. We find that the cause for this discrepancy is the method of assimilating the measured vertical shear into the simulation environment. The power law profile of Eq. (9) was fitted to measurements at and below hub height, and the shear coefficient was assumed to be independent of height. However, under stable conditions, the boundary layer height can be as low as 300 m. Thus, the vertical profile of wind speed between hub height (163 m) and rotor top (281 m) is not necessarily the same as between rotor bottom (45 m) and hub height. By analysing the vertical wind profile measured from a nacelle lidar, we found that the power law overestimates the mean wind speeds above hub height by 1–2 m s⁻¹. Consequently, the blade root bending moment had a larger range in the simulations and the damage equivalent moments were overestimated.

Figure 10 The damage equivalent moment at the root of the blade in the flapwise direction, as measured on the prototype (in black) and simulated (in red) in the aeroelastic solver by assimilating the measurements of inflow via spectral fitting with the Mann turbulence model, along with a power law profile for the vertical variation in the wind speed.

Table 6 The z scores for the damage equivalent moments on the root of the blade in the flapwise direction.

4 Discussion and conclusions

This study aimed to evaluate the accuracy of load simulations in which inflow turbulence is prescribed using the Mann model (Mann, 1994), by comparing the resulting fatigue loads with measurements from a 15 MW wind turbine prototype. Such measurements are rarely available or published, making this analysis a valuable contribution. In addition, the study provides insights into the influence of atmospheric stability on the fatigue damage accumulated by turbine towers and blades. The following conclusions can be drawn:

1. Analysis of the 1 Hz damage equivalent moment measured on the prototype demonstrated that atmospheric stability has a substantial effect on tower fatigue loads, while its impact on blade loads is less pronounced. These findings are consistent with Sathe et al. (2013), who additionally showed that the IEC recommendation can lead to estimations of lifetime fatigue damage that are up to 96 % higher than those obtained when accounting for the joint distribution of mean wind speed and stability at a given site. Our results therefore reinforce the case for including atmospheric stability in fatigue load assessments as a means of reducing over-design and associated costs. The alternative to the IEC approach is a site-specific fatigue load estimation, but current implementations typically assume neutral conditions to be dominant. While this assumption holds for many onshore sites, it may not be valid offshore or in coastal environments. Insufficient consideration of atmospheric stability at such sites could lead to the underestimation of lifetime fatigue damage on wind turbine towers.

2. A one-to-one comparison of measured and simulated fatigue loads, where the inflow was assimilated by fitting turbulence spectra to the Mann model, showed good overall agreement. The innovative aspect of this work is that simulations used the same controller as the prototype, ensuring a more realistic comparison. The damage equivalent moments from simulations were generally higher than measurements, but the margin of overestimation was small enough to be acceptable for engineering purposes. Improved inflow assimilation, for example through turbulence boxes constrained by nacelle lidar measurements, could further reduce these discrepancies. Alternatively, these differences could stem from uncertainties in the aerodynamic or structural properties of the simulated turbine, which are outside the scope of this study. Overall, we conclude that the Mann turbulence model is suitable for aeroelastic simulations of multi-megawatt wind turbines. Known limitations of the model, such as its assumption of Gaussian small-scale statistics and homogeneous turbulence, do not appear to significantly affect fatigue load predictions for tower and blade loads of a stand-alone turbine. Although the model was originally formulated for surface-layer, shear-driven turbulence, it successfully captures the key turbulence characteristics (length scale, turbulence intensity and coherence) most relevant for turbine loads. However, it is uncertain whether these results can be extrapolated to yaw and tilt moments on the main shaft.

3. In some one-to-one comparisons, significant overestimations of tower and blade fatigue loads were observed under certain stable conditions. For the tower, this was attributed to the Mann model overestimating vertical coherence in the along-wind component, while blade loads were overestimated due to an inaccurate shear profile. The physical mechanisms leading to the low measured coherence remain unclear, but future work could investigate whether the model of Chougule et al. (2018) is more accurate under some stable conditions. Furthermore, better representation of the shear profile above hub height in the simulations could be achieved through nacelle- or ground-based lidars capable of measuring up to the rotor top or by adopting an extended boundary layer wind profile (Gryning et al., 2007). Accurate assimilation of the shear profile can also be important for power curve verification (Wharton and Lundquist, 2012).

Subsequent research can be directed towards conducting one-to-one comparisons for wind turbines operating in wind farms. As the Mann model is assumed to describe atmospheric turbulence in time-domain simulations of wake meandering (Larsen et al., 2008), such studies would help to determine whether the model remains valid in this context or if extended formulations (Syed and Mann, 2024) are preferable. In addition, one-to-one comparisons could be performed for floating wind turbines, with inflow characterized using the Mann model, the Kaimal model (Kaimal et al., 1972) with exponential coherence (Davenport, 1961) or large-eddy simulations (Doubrawa et al., 2024).