Structural Reliability of Offshore Wind Turbines
Structural Reliability of Offshore Wind Turbines
Structural Reliability of Offshore Wind Turbines
Committee:
Robert B. Gilbert
Spyros A. Kinnas
Edward J. Powers
John L. Tassoulas
by
DISSERTATION
Presented to the Faculty of the Graduate School of
The University of Texas at Austin
in Partial Fulfillment
of the Requirements
for the Degree of
DOCTOR OF PHILOSOPHY
God,
the Reason of all reasons,
and
to my parents.
Acknowledgments
First and foremost, I sincerely thank my advisor, Dr. Lance Manuel, for
his guidance, encouragement, and unlimited availability for this research. I thank
members of my doctoral committee, Drs. John Tassoulas, Spyros Kinnas, Robert
Gilbert and Edward Powers, for their time and helpful comments. I also thank Dr.
Niels Jacob Tarp-Johansen for serving on my committee.
The financial support from National Science Foundation for this research is
gratefully acknowledged. Besides, I would like to acknowledge the financial support
provided by the Department of Civil, Architectural and Environmental engineering
during the first year of my doctoral program.
I thank Korn Saranyasoontorn, Jun Won Kang and Isamu Ikeda, Jr., my fellow
graduate students, for their help regarding solving several problems.
iv
Structural Reliability of Offshore Wind Turbines
Publication No.
v
bine model, we compare empirical short-term load distributions when two alternative
models for extremes—global and block maxima—are used. We develop a convergence
criterion, based on controlling the uncertainty in rare load fractiles, which serves to
assess whether or not an adequate number of simulations has been performed.
vi
Table of Contents
Acknowledgments iv
Abstract v
List of Tables xi
Chapter 1. Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Design Guidelines for Offshore Wind Turbines . . . . . . . . . . . . . 2
1.3 Statistical Load Extrapolation . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Challenges for Extrapolation of Wind Turbine Loads . . . . . . . . . . 5
1.4.1 Models for Extremes . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2 Robust Short-term Distributions . . . . . . . . . . . . . . . . . 7
1.4.3 Efficient Methods for Extrapolation . . . . . . . . . . . . . . . 7
1.4.4 Simulation Models . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Research Objectives and Methodology . . . . . . . . . . . . . . . . . . 8
1.5.1 Offshore Environment at Blyth and its Effects on Offshore Wind
Turbine Support Structure Loads . . . . . . . . . . . . . . . . . 9
1.5.2 Establishing Robust Short-term Load Distributions for Extrap-
olation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.3 Efficient Methods for Extrapolation of Long-term Offshore Wind
Turbine Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.4 Influence of Nonlinear Irregular Waves on Turbine Loads . . . . 11
1.6 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . 13
vii
2.3.2 Turbine Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 Representative Short-Term Response . . . . . . . . . . . . . . . . . . 33
2.5 Long-term Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5.1 Uncertainty in Long-term Load Estimates . . . . . . . . . . . . 40
2.5.2 Importance of Winds from the Shore . . . . . . . . . . . . . . . 46
2.6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 48
viii
4.4 Turbine Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.5 Short-term Distributions for Turbine Loads . . . . . . . . . . . . . . . 110
4.6 Long-term Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.6.1 The 2-D Environmental Contour (EC) Method . . . . . . . . . 114
4.6.2 Correction to the EC Long-term Loads . . . . . . . . . . . . . . 117
4.6.3 3-D Inverse FORM . . . . . . . . . . . . . . . . . . . . . . . . . 118
4.6.4 Control Actions and Number of Simulations . . . . . . . . . . . 121
4.6.5 Comparison of Alternative Models . . . . . . . . . . . . . . . . 124
4.6.5.1 POT and Global Maxima . . . . . . . . . . . . . . . . 125
4.6.5.2 Block Maxima . . . . . . . . . . . . . . . . . . . . . . . 127
4.6.5.3 Convergence Criterion for Required Number of Simula-
tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
4.7 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 130
ix
Chapter 6. Summary and Conclusions 193
6.1 Summary of Research Objectives . . . . . . . . . . . . . . . . . . . . . 193
6.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
6.2.1 Offshore Environment at Blyth and its Effect on Offshore Wind
Turbine Support Structure Loads . . . . . . . . . . . . . . . . . 194
6.2.2 Establishing Short-term Load Distributions for Extrapolation . 195
6.2.3 Efficient Methods for Extrapolation of Offshore Wind Turbine
Long-term Loads . . . . . . . . . . . . . . . . . . . . . . . . . . 197
6.2.4 Influence of Nonlinear Irregular Waves on Turbine Loads . . . . 198
6.3 Suggestions for Future Research . . . . . . . . . . . . . . . . . . . . . 200
Bibliography 202
Vita 216
x
List of Tables
2.1 Number of usable ten-minute summary data sets and ten-minute cam-
paign data sets from measurements at the Blyth site. . . . . . . . . . 20
2.2 Description of selected storms presented in Figs. 2.5 and 2.6. . . . . . 23
2.3 Distribution of datasets by wind-direction sectors. . . . . . . . . . . . 32
2.4 Statistics of the environment and load time series for the three selected
campaigns studied in Figs. 2.13, 2.14 and 2.15. . . . . . . . . . . . . . 34
3.1 Suggested block sizes (in seconds) for independent block maxima based
on mean (µB ) and mean plus one standard deviation (µB + σB ) values
from 200 simulations and tested at the 1% significance level. . . . . . 69
3.2 Average values of B computed using maxima from different block sizes
based on 20 simulations of four different loads for the most important
wind speed bins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
th
3.3 Estimates of the 84 percentile global maximum load for four different
load types as obtained from block maxima distributions using different
block sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.4 Estimates of the normalized 90% confidence interval on the 84th per-
centile load for OoPBM, FATBM, IPBM and OoPTD based on global
maxima from 30 simulations . . . . . . . . . . . . . . . . . . . . . . . 77
3.5 Estimates of the normalized 90% confidence intervals on the 84th per-
centile global maximum load for four different load types obtained from
different numbers of simulations. . . . . . . . . . . . . . . . . . . . . . 78
3.6 Estimates of the normalized 90% confidence intervals on the 84th per-
centile global maximum load for four different load types as obtained
from block maxima distributions using different block sizes. . . . . . . 80
xi
4.6 Effect of threshold level on the estimate of load fractiles for the design
environmental state for OoPBM (V = 14 m/s, Hs = 5.5 m, where
p3 = 0.99998997) and TBM (V = 16 m/s, Hs = 5.5 m, where p3 =
0.99999613). Threshold = Mean + Nσ ×(Standard deviation). . . . . 126
4.7 Average plus one standard deviation values of B from block maxima
with different block sizes based on 60 simulations for the out-of-plane
bending moment at the blade root (OoPBM) and based on 150 simu-
lations for the fore-aft tower bending moment at mudline (TBM). . . 127
4.8 Estimates of normalized 90% confidence interval on the 84th percentile
load for design environmental states for the out-of-plane bending mo-
ment at the blade root (OoPBM) and the fore-aft tower bending mo-
ment at the mudline (TBM). . . . . . . . . . . . . . . . . . . . . . . . 129
5.1 Statistical moments of the linear and nonlinear sea surface elevation
processes simulated for a JONSWAP spectrum with Hs = 12 m and
Tp = 14 sec. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5.2 Comparison of statistics of the bending moment at the base of the rigid
stand-alone monopile, averaged over fifty simulations, when linear and
nonlinear waves are simulated using a JONSWAP spectrum with Hs
= 5.5 m and Tp = 11.2 sec . . . . . . . . . . . . . . . . . . . . . . . . 153
5.3 Comparison of hydrodynamic loads on a stand-alone monopile ana-
lyzed quasi-statically and water particle kinematics, averaged over fifty
simulations, computed with linear and nonlinear irregular waves using
a JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec. . . . . . 156
5.4 Comparison of statistics, averaged over fifty simulations, of bending
moment at the base of the rigid stand-alone monopile due to drag and
inertia forces, when linear and nonlinear waves are simulated using a
JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec. . . . . . . 158
5.5 Comparison of statistics of Morison forces (per unit length) at the mean
sea level, averaged over fifty simulations, with statistics obtained from
the theoretical model, when linear irregular waves are simulated with
a JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec. . . . . . 161
5.6 Values of steepness, variance and skewness of the nonlinear sea surface
elevation process for various choices of spectral peak period, Tp , for a
significant wave height, Hs , of 5.5 m, when waves are simulated using
a JONSWAP spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . 162
5.7 Ratio of the maximum quasi-static base moment on a stand-alone
monopile computed using nonlinear waves to that using linear waves
(based on 50 ten-minute simulations). . . . . . . . . . . . . . . . . . . 168
5.8 Ten-minute maximum fore-aft tower bending moment (in MN-m) at
mudline, averaged over 20 simulations, computed with linear and non-
linear irregular waves using a JONSWAP spectrum with Hs = 5.5 m
and Tp = 11.2 sec, and V = 16 m/s. . . . . . . . . . . . . . . . . . . 170
5.9 Comparison of statistics of the sea surface elevation process simulated
using linear and nonlinear irregular wave models for a JONSWAP spec-
trum with Hs = 7.5 m and Tp = 12.3 sec. . . . . . . . . . . . . . . . . 172
xii
5.10 Comparison of statistics, averaged over 50 simulations, of the fore-aft
tower bending moment at the mudline for linear and nonlinear waves,
for a JONSWAP spectrum with Hs = 7.5 m and Tp = 12.3 sec, and
wind speed, V = 16 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . 173
5.11 Largest observed ten-minute maximum fore-aft tower bending moment
at the mudline, from 500 simulations, for the governing environmental
state for the linear wave model (V = 16 m/s, Hs = 5.5 m) and for the
nonlinear wave model (V = 18 m/s, Hs = 7.5 m). . . . . . . . . . . . 181
5.12 Estimates of the 84th percentile 10-min maximum fore-aft tower bend-
ing moment at the mudline from 10 separate sets of 50 simulations
each. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
5.13 Estimates of the 90% relative confidence bound (RCB) on the 84th per-
centile 10-min maximum fore-aft tower bending moment at the mudline.184
5.14 Extrapolated loads (10-min maximum fore-aft tower bending moment
at mudline) from 50 simulations. . . . . . . . . . . . . . . . . . . . . . 186
xiii
List of Figures
xiv
2.13 Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the
mudline bending moment for the campaign time series, ShoreV18H0N1 36
2.14 Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the
mudline bending moment for the campaign time series, ShoreV18H0N9 36
2.15 Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the
mudline bending moment for the campaign time series, SeaV12H3N41 38
2.16 Comparison of probability of load exceedance curves for winds from
the sea, winds from the shore, and winds from all directions. . . . . . 39
2.17 Long-term loads (mudline bending moments) for 1-year and 20-year
return periods for different wind-direction sectors. Also shown are 90%-
confidence intervals along with mean values for the ‘with-uncertainty’
case based on bootstrap methods applied to the limited data in each
sector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.18 Comparison of 90-percent confidence intervals and mean value esti-
mates of the 20-year load based on different methods of estimation of
parameters of the short-term load distribution: (a) regression over the
above-median loads data; (b) regression over all the loads data; (c) col-
location using the 50th and 90th load percentiles; and (d) the method
of moments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.19 Comparison of 20-year loads when the environmental random variable
vector for each wind-direction sector, X, consists of (a) mean wind
speed (V ) and significant wave height (Hs ); and (b) mean wind speed
(V ) and turbulence standard deviation (σV ). . . . . . . . . . . . . . . 48
xv
3.7 Short-term distributions for OoPBM, FATBM, IPBM and OoPTD
based on 30 simulations. Bins selected have largest 90% relative confi-
dence bounds on the 84th percentile load as computed based on boot-
strap techniques. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
xvi
4.14 Load distributions of the POT data for governing environmental states
based on (a) 60 simulations of OoPBM for a mean wind speed of 14
m/s and a significant wave height of 5.5 m; and (b) 150 simulations of
TBM for a mean wind speed of 16 m/s and a significant wave height
of 5.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.15 Two-parameter Weibull distribution fits to POT data for (a) OoPBM
(for V = 14 m/s, Hs = 5.5 m), and based on 60 simulations (b) TBM
(for V = 16 m/s, Hs = 5.5 m) and based on 150 simulations. . . . . . 124
4.16 Two-parameter Weibull distribution fits to global maxima data for (a)
OoPBM (for V = 14 m/s, Hs = 5.5 m), and based on 60 simulations
(b) TBM (for V = 16 m/s, Hs = 5.5 m) and based on 150 simulations. 125
5.1 Simulation of sea surface elevation using a delta function power spec-
trum showing differences between linear and nonlinear waves. . . . . . 148
5.2 Time series of the sea surface elevation process simulated using a JON-
SWAP spectrum with Hs = 12 m and Tp = 14 sec. . . . . . . . . . . 149
5.3 Probability of exceedance of different sea surface elevations simulated
using linear and nonlinear wave theories for a JONSWAP with Hs =
12 m and Tp = 14 sec. The data are from 100 simulations and the
Gaussian distribution is based on the computed mean and standard
deviation with the linear wave model. . . . . . . . . . . . . . . . . . . 151
5.4 Power spectrum of (a) sea surface elevation and (b) bending moment at
the base of the rigid stand-alone monopile, when linear and nonlinear
waves are simulated using a JONSWAP spectrum with Hs = 5.5 m
and Tp = 11.2 sec. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
5.5 Variation of the ten-minute maximum bending moment at the base
of a rigid stand-alone monopile, averaged over fifty simulations, with
wave steepness, when linear and nonlinear waves are simulated using
a JONSWAP spectrum with Hs = 5.5 m while Tp values are chosen
according to the wave steepness. . . . . . . . . . . . . . . . . . . . . . 163
5.6 Variation of the ten-minute maximum horizontal water particle accel-
eration at the mean sea level, averaged over fifty simulations, with
wave steepness, when linear and nonlinear waves are simulated using
a JONSWAP spectrum with Hs = 5.5 m while Tp values are chosen
according to the wave steepness. . . . . . . . . . . . . . . . . . . . . . 164
5.7 Variation of the second-order transfer function, D+ , for sum-frequency
interactions, (ω + ω) for a water depth of 20 meters. . . . . . . . . . . 167
5.8 (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for wind speed (for V = 16 m/s) and
of sea surface elevation (for linear and nonlinear waves simulated using
a JONSWAP spectrum with Hs = 7.5 m and Tp = 12.3 sec). . . . . . 176
5.9 (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for fore-aft tower base shear (FATBS)
and tower bending moment (FATBM) at the mudline (for Hs = 7.5 m
and Tp = 12.3 sec, when wind is not included). . . . . . . . . . . . . . 177
xvii
5.10 (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for fore-aft tower base shear (FATBS)
and tower bending moment (FATBM) at the mudline (for Hs = 7.5 m
and Tp = 12.3 sec, when wind is included, i.e., V = 16 m/s). . . . . . 178
5.11 Empirical probability distributions of ten-minute maxima of the fore-
aft tower bending moment at the mudline based on 500 ten-minute
simulations for (a) V = 16 m/s, Hs = 5.5 m; and (b) V = 18 m/s and
Hs = 7.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
5.12 Empirical probability distributions of ten-minute maxima of the wave
height based on 500 ten-minute simulations for (a) V = 16 m/s, Hs =
5.5 m; and (b) V = 18 m/s and Hs = 7.5 m. . . . . . . . . . . . . . . 182
5.13 Ten-minute time series of the wave elevation, wind speed, and fore-aft
tower bending moment showing almost concurrent occurrences of the
maximum wave height and the maximum tower load for V = 18 m/s,
Hs = 7.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
5.14 Scatter plots showing whether the occurrences of the maximum fore-
aft tower bending moment in ten minutes and the occurrence of the
maximum instantaneous wave height are close (i.e., well correlated) or
far apart (i.e., less correlated). Plots are for V = 16 m/s and Hs = 5.5
m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
5.15 Scatter plots showing whether the occurrences of the maximum fore-
aft tower bending moment in ten minutes and the occurrence of the
maximum instantaneous wave height are close (i.e., well correlated) or
far apart (i.e., less correlated). Plots are for V = 18 m/s and Hs = 7.5
m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
xviii
Chapter 1
Introduction
1.1 Background
Wind energy, growing annually at a rate of 20-30% [66] over the last decade, is
the fastest growing energy resource globally. The main drivers for such growth are an
increased concern for global warming, an increase in the prices of oil and gas, and a
decrease in the cost of wind energy due to technological developments. Wind energy
accounts for about 3% of the total installed electricity capacity in Europe, with a
target of about 12% by the year 2020, of which about one-third, or about 30 GW
(gigawatts), is expected to come from offshore projects [21, 22]. While all the installed
wind power capacity (about 17 GW by 2007 [1]) in the United States is from onshore
developments, the potential for offshore wind power in the U.S., in water depths less
than 30 meters where current fixed-bottom support structure technologies can be
used, is about 98 GW [70]. Clearly, offshore wind energy is becoming an important
part of the overall energy mix outside the U.S. and has great potential within the
U.S. as well.
1
required for their design, becomes more difficult, mainly due to consideration of wave
forces in addition to those from the wind. In fact, even 17 years after the first offshore
wind farm was commissioned, the design guidelines for the offshore wind turbines are
still a state of in development, and there is a clear need for additional research in this
area.
2
loads for offshore wind turbines is our primary interest in this dissertation.
Design Load Case (DLC) 1.1b of the IEC 61400-3 draft design guidelines [43]
for offshore wind turbines, which is based on DLC 1.1 of the IEC 61400-1 guide-
lines [42] for onshore wind turbines, recommends the use of statistical extrapolation
to predict turbine nominal loads for an ultimate limit state. Direct integration is
the primary method for extrapolation, in which one estimates the turbine load, lT ,
associated with a target probability of exceedance, PT , or, equivalently, with a target
return period of T years, as follows:
Z
PT = P [L > lT ] = P [L > lT |X = x] fX (x)dx (1.1)
X
where fX (x) represents the joint probability density function of the environmental
random variables, X, and L is the random variable representing the load of interest.
For different trial values of the load, lT , Eq. 1.1 enables one to compute the long-term
probability by integrating the short-term load exceedance probability conditional on
X, i.e., P [L > lT |X = x], with the relative likelihood of different values of X. The
load level at which the computed long-term probability matches the target proba-
bility is the desired nominal load. The direct integration method, while exact, is
computationally expensive as one is required to integrate over the entire domain of
all the environmental random variables. For offshore wind turbines, X is assumed
to be comprised primarily of the ten-minute average wind speed, V , at hub height
in the along-wind direction and the significant wave height (four times the standard
deviation of the sea surface elevation process), Hs , for waves assumed to be aligned
with the wind.
3
Given
X=x
r=r+1
no All x
realized?
Figure 1.1: Flow chart representing the steps involved in establishing short-term and
long-term load distributions based on turbine response simulations.
The flow chart in Fig. 1.1 shows how turbine response simulations are used to
establish the long-term load distribution using the direct integration approach. For
a given value of environmental random variables, X = x, we first seek to establish
the short-term distribution of load conditional on the environment, FL|X =x (l). To
do so, M load extremes, Lr (r = 1, ..., M ; where M is usually a number not smaller
than 6), are required. These can then be used to establish an empirical (short-
term) distribution function, to which a parametric probability distribution model
may be fitted if desired. A global maximum, which is the maximum value from a
ten-minute (the chosen duration here, as recommended by design standards [42, 43])
time series of turbine response, is an example of the load extreme. To obtain one
time series realization of the stochastic turbine response process, we first simulate a
4
ten-minute realization of the inflow stochastic wind field on the rotor plane, which
is possible using turbulence simulation software such as TurbSim [45]. The incident
wave process on the support structure is also simulated, commonly using a linear
irregular (stochastic) wave modeling approach. Once these environmental inputs are
in place, the response of the turbine is computed using an aero-servo-hydro-elastic
analysis in the time domain. This is usually done using available turbine response
simulation software such as FAST [47]. By changing the set of random seeds, sr , that
produces a different realization of the stochastic wind and wave processes for the same
mean environmental conditions x, and then generating a different realization of the
turbine response, one can obtain the desired M load extremes. Once the short-term
load distributions for all X (i.e., for all possible combinations of mean wind speed
and significant wave height) are established, we can integrate them with the joint
probability distribution function of the environmental random variables, fX (x), to
obtain the long-term load distribution using Eq. 1.1.
5
to active control systems that adjust the pitch angle of the turbine blades depending
on the inflow wind speed, the turbine system effectively changes in real time. As a
result, the structural response of a turbine can show significant variability between
two random simulations even for the same mean environmental state. In addition, the
extrapolation of turbine loads needs to recognize the dependence on two (or more)
random processes representing the environment — wind and waves, say — each of
which can influence turbine loads in distinct ways. In fact, each box in the flow chart
shown in Fig. 1.1 poses separate questions for research, some of which are discussed
below.
The extremes from a time series of any turbine load for a given duration can be
extracted in several ways [15]. The global maximum, which is the single largest value
from a time series, is the most common and the simplest method, and is discussed
in detail in the literature [30, 55]. Another extreme model is the peak-over-threshold
(POT) method, in which the maximum value from each segment of a time series that
lies between two successive upcrossings of a chosen threshold is retained as a load
extreme. Yet another model for extremes is to use a block maximum, in which one
partitions the time series into individual non-overlapping blocks of constant duration,
and the largest values from each of these blocks constitute a set of block maxima.
There are several important questions regarding extrapolation using different the
extreme models, which have not yet been addressed in the wind turbine design guide-
lines [42, 43]. These questions include: whether the long-term loads predicted using
different extremes models are comparable; whether the additional extreme statistics
from methods like POT and block maxima (over that from global maxima) can re-
duce the need for large number of simulations; and, lastly, how to choose an optimal
6
threshold or an optimal block size, if one exists, for the POT and block maxima
methods, respectively. Of these questions, only those regarding the POT methods
have been addressed in recent research [83].
The direct integration method for statistical extrapolation requires that re-
peated computationally expensive simulations be run over the entire domain of en-
vironmental random variables. The simulation effort, even with a limited number
of simulations, can involve hundreds of hours of computer time on today’s personal
computers. It is thus of interest to explore alternative extrapolation techniques that
are more efficient, such as the inverse first-order reliability method [107], which has
been demonstrated to be accurate for onshore wind turbines [89], for possible use in
predicting long-term loads on an offshore wind turbine.
7
1.4.4 Simulation Models
We have four main research objectives, all related to the prediction of long-
term loads for offshore wind turbines using statistical extrapolation: (1) to understand
the nature of the offshore environment at a typical offshore wind farm site and how
the wind-wave environment there can influence long-term loads on the wind turbine
support structure; (2) to develop an improved understanding of extrapolation proce-
dures for wind turbine extreme loads and suggest guidelines for their practical use;
(3) to assess the accuracy of the inverse first-order reliability method, an efficient but
8
approximate method, to predict long-term loads; and (4) to investigate the effect of
a second-order nonlinear irregular wave model on wind turbine extreme loads. In the
following, some details on how we intend to achieve these objectives are presented.
We broadly classify winds at the Blyth site based on whether they are blow-
ing from the shore towards the sea or from the sea towards the shore; we further
subdivide these two wind regimes into 30◦ wind-direction sectors and examine the
wind-wave correlation more closely and how the loads are influenced by wind and
waves differently in each sector. In particular, we establish (conditional) short-term
distributions of the load extremes (ten-minute global maxima) given wind and waves
for each sector, and then derive long-term loads by accounting for the likelihood of
different wave heights and wind speeds in each sector. We also compute confidence
bounds on long-term load predictions using the bootstrap technique.
9
1.5.2 Establishing Robust Short-term Load Distributions for Extrapola-
tion
10
threshold method is used in conjunction with inverse FORM.
11
predictions that arise from the use of the alternative wave models.
1.6 Limitations
We focus only on ultimate limit states, while ignoring fatigue loads that can
sometimes drive the overall design of some wind turbine components. Moreover, we
consider loads only during the operating state of a wind turbine while severe loads may
also result during parked, start-up, and shut-down states, which we ignore. All our
simulations, and the resulting conclusions, are based on one offshore wind turbine
model—the NREL 5MW baseline offshore wind turbine. While this model closely
represents modern utility-scale offshore wind turbines, our conclusions may not be
12
valid for all wind turbines. Regarding the environment, we assume that the wind
models described in the design guidelines, which are based on models for onshore
sites, are adequate for offshore sites as well; this may not strictly be true. As far as
waves are concerned, we investigate the influence of a second-order nonlinear irregular
wave model on turbine loads, compared to the commonly used but less realistic linear
irregular wave model. However, other wave models may be needed to describe higher
order nonlinearities. Besides, we do not address breaking of waves (which is generally
thought to be important) since our focus is on long-term probabilistic load prediction
by simulation and there is no well-established way to model irregular breaking waves.
Lastly, in line with current practice, we assume that the monopile foundation is
rigidly connected at the base, which is not likely to be the case unless the seabed at
an offshore wind turbine site consists of rock-like conditions.
One should interpret the conclusions from this dissertation in light of the
above limitations. Future research may be needed to address the influence of these
limitations on predicted long-term loads for offshore wind turbines.
The organization of this dissertation follows the outline of the research objec-
tives and methodology. In Chapter 2, we study characteristics of the wind and waves
at the Blyth site, and investigate their influence on long-term loads on the support
structure of a 2MW offshore wind turbine at the Blyth wind farm. The remaining
three chapters concern simulation studies on the NREL 5MW wind turbine model.
In Chapter 3, we use a data set of simulated loads on the 5MW onshore turbine
to develop procedures for establishing robust short-term distributions using global
maxima and block maxima methods, while also investigating the independence of
13
selected block maxima. In Chapter 4, we assess the accuracy of the inverse first-order
reliability method, a more efficient method for extrapolation than direct integration,
in predicting long-term loads for the 5 MW offshore wind turbine, while highlighting
the importance of variability in turbine loads. We also compare long-term loads pre-
dicted using the three extreme models — global maxima, POT, and block maxima. In
Chapter 5, we introduce the second-order nonlinear irregular wave model, we present
a procedure for efficient simulation of such waves, and we compare extreme loads on
the turbine support structure due to linear and nonlinear irregular waves. We also
discuss in detail characteristics of the wave loading due to alternative wave models.
Finally, in Chapter 6, we summarize conclusions from our research.
14
Chapter 2
2.1 Introduction
15
Figure 2.1: Turbines at the Blyth site.
water site about one kilometer off the Northumberland shore in northeast England. At
this site, a 2MW wind turbine was instrumented, for which data on the environment
and loads were recorded over a period of roughly sixteen months [10]. From this
measurement campaign, loads data available are in two different formats: as ten-
minute statistics (referred to as “summary” data) and as full time series (referred
to as “campaign” data). An interesting feature of this site is that there are starkly
different environmental characteristics associated with winds blowing from various
directions—for instance, in the degree of correlation between wind speed and wave
height and in the turbulence characteristics of the wind. As a result, the turbine
response is expected to be different for the different wind-direction sectors. The winds
16
at this site are broadly classified based on whether they are blowing from the shore
towards the sea or from the sea towards the shore. We further subdivide these two
wind regimes into 30◦ wind-direction sectors and examine more closely how the loads
are influenced by wind and waves differently in each sector. In particular, we seek
to establish (conditional) short-term distributions of the load given wind and waves
for each sector, and then derive long-term loads by accounting for the likelihood of
different wave heights and wind speeds in each sector. We use the summary data,
which are ten-minute global maxima, to establish short-term and long-term load
distributions. We also use the campaign data to study characteristics of the inflow
environment and the turbine response in some detail using time series and power
spectra.
Lessons learned from this study of field data on wind turbine loads especially
regarding nonlinear characteristics of waves and the need for efficient extrapolation
procedures in long-term load predictions serve to motivate the work in rest of this
dissertation.
The Blyth project is an experimental wind farm consisting of two 2MW Vestas
V66 wind turbines. The Blyth site is located on the northeastern coast of England, off
the Northumberland shore. The wind turbines are located approximately 1 km from
the shoreline. The mean water depth at the instrumented turbine varies between
a Lowest Astronomical Tide (LAT) level of 6 m and a Mean High Water Springs
(MHWS) level of 11 m. The average water depth at the turbine location is approx-
imately 9 m. One of the two turbines (the southern turbine in Fig. 2.2a) at Blyth
was instrumented as part of a research project funded by the European Commission;
17
(a) (b)
Figure 2.2: (a) Location of the turbines and an onshore meteorological mast at the
Blyth site; (b) Layout of the turbine instrumentation (from Camp et al. [10])
it has a hub height of 62 m above the LAT level and a blade diameter of 66 m. The
turbines are sited on sharply sloping submerged rock, known as the ‘North Spit,’ in
rock-socket type foundations. This local bathymetry results in rather large waves,
including breaking waves, at the turbine.
Field measurements were collected for sixteen months between October 2001
and January 2003; thus, the period of data collection covered more than one full
winter season. Measured data included wind speed and direction at the nacelle,
sea surface elevation, and bending moments at several vertical stations along the
tower and the pile. One of these stations—the mudline bending moment—is our load
variable of interest here. (Blade loads were also recorded but are not analyzed in the
present study.) Figure 2.2b shows the layout of the turbine instrumentation. The
hub-height wind speed was measured using an anemometer placed on the turbine
nacelle, downwind of the rotor, and thus the wind speed data were influenced by
the turbine wake. The hub-height wind speed data were calibrated [10] such that
the mean wind speed measured at the nacelle was approximately equal to the mean
18
free wind speed measured at a similar elevation at a nearby shore location. No
corrections were applied to the nacelle’s turbulence measurements, and appropriate
turbulence characteristics may only be obtained from the wind speed data measured
at a meteorological mast located on the shore, approximately 1 km from the turbine.
Sea surface elevation data were measured using a wave radar system located at the
entrance platform of the turbine, 11.7 m above the LAT line. Additional details
regarding the data and measurement system may be found in the report by Camp et
al [10].
Time series data in ten-minute segments were sampled at 40 Hz, and mini-
mum, maximum, mean, and standard deviation values for each channel were recorded
as part of the statistics comprising the ‘summary’ data sets. In addition, time his-
tories of ten-minute duration were also recorded when a predetermined set of trigger
conditions were met; these time-series data are referred to as the ‘campaign’ data.
The triggers for normal operating conditions were related to variations in mean wind
speed and direction, significant wave height, and tidal water level. Several other trig-
ger conditions were defined in terms of the turbine state (operating, idling, start-up,
and shut-down) and extreme environmental events that included wind speeds above
25 m/s and significant wave heights above 4 m [10].
In this study, we will rely primarily on the statistics from the summary data to
derive distributions for the environmental random variables and the mudline bending
moment. While almost 64,000 data sets or ten-minute samples were continuously
recorded as part of the summary data, only about 2,300 sets were found to be usable
after data sets with bad or missing values were discarded. Similarly, out of almost
2,900 ten-minute time series recorded in the campaign data, only 1,014 were found
to be suitable for providing meaningful statistics. It should be noted that the criteria
19
Table 2.1: Number of usable ten-minute summary data
sets and ten-minute campaign data sets from measure-
ments at the Blyth site.
Winds from Winds from Winds from
Data Type
the sea the shore all directions
Summary 1,398 903 2,301
Campaign 646 368 1,014
for “good” datasets included the requirement for: (1) simultaneously available values
for statistics/parameters of interest (mean wind speed, turbulence standard devia-
tion, significant wave height, and mudline bending moment); and (2) realistic values
of physical variables (e.g., data with zero or negative wind speeds were omitted).
Furthermore, we focussed only on datasets recorded while the turbine was operating.
The distribution of the usable data of the two types is summarized in Table 2.1.
2.3.1 Environment
20
random variable vector, X, has components (Θ, V, Hs ). To establish the joint distri-
bution of the environment, we will estimate the marginal distribution of mean wind
direction, the distribution of mean wind speed conditional on mean wind direction,
and the distribution of significant wave height conditional on both the mean wind
direction and mean wind speed.
We first discuss the distribution of mean wind direction. The wind direction
here refers to the direction (angle in degrees) from which winds are blowing. We
divide the summary data into 30◦ wind-direction sectors. The wind direction at the
site ( Fig. 2.2a) is measured clockwise from geographical north, which corresponds to
a wind direction of 0◦ . It is useful to classify the data according to wind directions
into two broad regimes [10]—namely, “winds from sea” and “winds from shore.” At
this site, wind directions between 0◦ and 150◦ constitute “winds from sea,” while
wind directions between 180◦ and 330◦ constitute “winds from shore.” The set of all
wind directions will be referred to as “all-direction” winds.
Number of Occurrences
200
150
100
50
0
0
60
120 24
θ( 180 20 s)
de 240 16 /
gr
ee 300360 12 , V (m
8 ed
)s 4 d Spe
0 Win
an
Me
(a) (b)
Figure 2.3: (a) Histogram showing number of occurrences or data sets for different
mean wind directions while the turbine is operating; (b) histogram showing division
of data according to both mean wind direction and mean wind speed.
21
directions while the turbine is operating; it is observed that winds from the sea, partic-
ularly those with wind directions between 60◦ and 120◦ , occur more often than winds
from the shore. Figure 2.3b shows a histogram summarizing the joint occurrences of
mean wind direction and mean wind speed. It is observed that winds from the sea are
associated with a greater relative frequency of higher wind speeds than lower wind
speeds. This is in contrast with winds from the shore for which lower wind speeds
are more common. This is also reflected in the respective mean wind speeds for each
data set separately; the mean wind speeds are 9.57 m/s, 6.25 m/s, and 7.70 m/s,
respectively, for winds from the sea, winds from the shore, and all-direction winds.
Figure 2.3a clearly shows that some wind-direction sectors are scarcely sampled sug-
gesting that some wind directions occurred rarely; furthermore, Fig. 2.3b shows that
many wind speed bins for a given wind direction have a small amount of data or, in
some cases, no data at all. This underlines the limited nature of the data available
for this study.
Before establishing probability distributions for the mean wind speed and sig-
nificant wave height, we study the correlation between wind and waves as seen in the
data. Figure 2.4 shows observed values of mean wind speed versus significant wave
height and suggests significant scatter in the data, with winds from sea and shore fol-
lowing different trends. The winds from the sea suggest greater correlation between
wind speed and wave height, which may be explained by the longer fetches typically
associated with winds from the sea. In the present study, we are only interested in
evaluating loads associated with the turbine in an operating state; as such we are not
interested in wind speeds below the cut-in wind speed of 4 m/s.
An interesting feature exhibited by the data for winds from the shore is the
presence of some large wave heights associated with small wind speeds below 4-5 m/s,
22
Winds from Sea
5
0
0 5 10 15 20 25
Mean Wind Speed, V (m/s)
Figure 2.4: Scatter diagram showing mean wind speed versus significant wave height.
indicating weaker correlation between wind and waves. The difference in correlation
in wind-wave data for winds from the shore versus winds from the sea may be partly
explained by the contrasting nature of storm profiles for the two wind regimes. To
better understand this, we select two storms—one associated with winds from the sea
and the other with winds from the shore; some key parameters that describe these
storms are summarized in Table 2.2. The durations of the storms for winds from the
sea and winds from the shore are three days and two days, respectively.
Table 2.2: Description of selected storms presented in Figs. 2.5 and 2.6.
Mean wind direction Storm dates No. of Max V Max Hs
(from) (degrees) from to samples (m/s) (m)
23
20 5 5
4 4
15
3 3
10
V→ Hs→
2 2
5
1 1
0 0 0
00:00 12:00 00:00 12:00 00:00 12:00 00:00 0 5 10 15 20 25
Time of day (hours) Mean Wind Speed, V (m/s)
(a) (b)
Figure 2.5: (a) Storm profile and (b) wind-wave scatter diagram for winds from the
sea. The duration of the storm shown is from Dec 7 to Dec 10, 2002.
20 5 5
V→ 4 4
15
Hs→
3 3
10
2 2
5
1 1
0 0 0
00:00 12:00 00:00 12:00 00:00 0 5 10 15 20 25
Time of day (hours) Mean Wind Speed, V (m/s)
(a) (b)
Figure 2.6: (a) Storm profile and (b) wind-wave scatter diagram for winds from the
shore. The duration of the storm shown is from Jan 22 to Jan 24, 2003.
Figure 2.5a shows the development of the selected storm associated with winds
from the sea; it is seen that both the mean wind speed and significant wave height
increase or decrease more or less concurrently. As a result, relatively strong positive
correlation is observed between the mean wind speed and significant wave height for
this storm as can be confirmed from Fig. 2.5b. On the other hand, the development
24
30 < θ < 60 60 < θ < 90 90 < θ < 120
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
Significant Wave Height, H (m)
0 0 0
0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25
s
120 < θ < 150 180 < θ < 210 210 < θ < 240
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
0 0 0
0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25
240 < θ < 270 270 < θ < 300 300 < θ < 330
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
0 0 0
0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25
Mean Wind Speed, V (m/s)
Figure 2.7: Wind-wave scatter diagrams for each of the 30◦ wind-direction, θ, sectors.
of significant wave heights generally lags the development of mean wind speeds for
the storm associated with winds from the shore, as shown in Fig. 2.6a. Consequently,
large significant wave heights are often observed with small mean wind speeds and vice
versa, resulting in relatively weaker correlation between the two, as may be observed
from Fig. 2.6b.
While winds from the sea and shore display different trends (in relation to the
accompanying waves), there is still significant scatter in each of these regimes sepa-
rately, as is clear from Fig. 2.4 and, therefore, we study the correlation between wind
and waves when the data are further divided into 30◦ wind-direction sectors. Fig-
ure 2.7 shows that the correlation between wind and waves depends on the direction
25
of the wind. This correlation can be physically understood to depend on whether,
in a given sector, the winds blow from the open sea, from and over land, or along
the shoreline. The wind-direction sectors between 30◦ and 120◦ correspond to winds
blowing from the open sea and, therefore, mean wind speeds and significant wave
heights for these directions are highly correlated due to the longer fetch associated
with such conditions. Winds from the sectors between 180◦ and 270◦ come from over
land passing over approximately 1 km of water before arriving at the turbine site;
associated waves are, thus, not expected to develop much due to the limited fetch.
This results in small wave heights for a given wind speed, for these wind directions.
The wind-direction sectors between 120◦ and 150◦ as well as those between 270◦ and
330◦ correspond to winds blowing roughly parallel to the shoreline and, thus, the en-
vironment for these sectors is influenced by characteristics of the land as well as the
sea. As a result, wave heights for the 120◦ -150◦ wind-direction sector are relatively
smaller among the entire set of cases for winds from the sea (i.e., in the 0◦ -150◦ range),
while wave heights for directions between 270◦ and 330◦ (and especially in the 300◦ -
330◦ sector) are relatively larger among the entire set of cases for winds from the shore
(i.e., in the 180◦ -330◦ range). These observations serve to illustrate that one should
recognize the different degrees of correlation between wind speed and wave height as
a function of the dominant wind direction as this correlation can have a bearing on
turbine loads when wind and wave loads are of comparable importance.
We assume that the ten-minute mean wind speed, V , conditional on the mean
wind direction, follows a Rayleigh distribution, which is often used to represent the
ten-minute mean wind speed’s probability distribution (see, for example, IEC design
guidelines [42, 43]). This Rayleigh distribution is truncated both below a cut-in
wind speed, Vin , of 4 m/s, and above a cut-out wind speed, Vout , of 25 m/s, since
26
Figure 2.8: Sample averages and sample 90%-confidence intervals (5- and 95-
percentile levels) of the ten-minute mean wind speed for different mean wind di-
rections.
we are interested only in operating turbine states. The expression for the truncated
cumulative distribution function (CDF), FV |Θ (v), of V is thus:
G(Vin ) − G(v)
FV |Θ (v) = (2.1)
G(Vin ) − G(Vout )
" µ ¶#2
v 2
G(v) = exp − ; α(θ) = √ µV |Θ (θ)
α(θ) π
where α is the scale parameter of the Rayleigh distribution, which is estimated from
the expected value, µV |Θ (θ), of V for a given mean wind direction, Θ.
Figure 2.8 shows the expected value (actually, the sample average) of the ten-
minute mean wind speed for each wind-direction sector along with sample 5- and
95-percentile non-exceedance probability levels. It is again seen that winds from the
sea (in the 0◦ -150◦ wind-direction sectors) are associated with somewhat higher mean
wind speeds than winds from the shore.
27
To represent variability in the wave climate, we assume that the significant
wave height, Hs , conditional on the mean wind speed, V , and the mean wind direction,
Θ, follows a Weibull distribution (as recommended in the DNV offshore standard [16]).
The expression for the cumulative distribution function of Hs conditional on V and
Θ, namely FHs |V,Θ (h), is given by:
" µ ¶k(v(θ)) #
h
FHs |V,Θ (h) = 1 − exp − (2.2)
λ(v(θ))
The shape parameter, k, and the scale parameter, λ, of the Weibull distribution
depend on the mean wind speed, and implicitly on the mean wind direction. The
data for each wind-direction sector are binned according to a moving window along
wind speed. The size of this wind speed window is 2 m/s; the window is moved in
increments of 0.5 m/s over the range from 4 m/s to 24 m/s. For each such window,
the Weibull parameters are estimated by the method of moments, using the following
equations:
v ¡ ¢
u
u Γ 1 + k2 λ
COVHs |V,Θ = t£ ¡ ¢¤2 − 1; µHs |V,Θ = ¡ ¢ (2.3)
Γ 1 + k1 Γ 1 + k1
where µHs |V,Θ and COVHs |V,Θ are the expected value and the coefficient of variation,
respectively, of significant wave height conditional on mean wind speed and mean wind
direction. As an illustration, estimates of these parameters are shown in Fig. 2.9 for
one wind-direction sector (90◦ -120◦ ) for winds from the sea and one wind-direction
sector (240◦ -270◦ ) for winds from the shore. It is observed that the shape parameter,
k, follows similar trends (increasing with wind speed) for both the wind-direction
sectors selected. This is because the coefficient of variation on significant wave height,
which is roughly inversely proportional to the shape parameter, is generally found to
28
Weibull Shape Parameter, k (m)
3
10
2
5
1
0 0
4 6 8 10 12 14 16 18 20 22 24 4 6 8 10 12 14 16 18 20 22 24
Mean Wind Speed, V (m/s) Mean Wind Speed, V (m/s)
(a) (b)
Figure 2.9: Weibull distribution parameters for significant wave height conditional
on mean wind speed and mean wind direction—(a) shape parameter, k; (b) scale
parameter, λ. Estimates are for the 90◦ -120◦ wind-direction sector representative of
winds from the sea and the 240◦ -270◦ wind-direction sector representative of winds
from the shore.
decrease with increasing wind speeds for both winds from the sea as well as from
the shore. On the other hand, the scale parameter, λ, is considerably larger for the
selected wind-direction sector for winds from the sea compared to that for the selected
wind-direction sector for winds from the shore. This is due to the fact that, for a
given wind speed, the expected value of the significant wave height, which is directly
proportional to the Weibull scale parameter, is larger for winds from the sea than for
winds from the shore.
We are interested in the mudline bending moment as our load variable. The
summary data provide the maximum (over ten minutes) mudline bending moment
recorded in two orthogonal directions, x and y, which are fixed in space. The square
root of the sums of squares (SRSS) of these two maximum values is always larger than
the maximum of the time-varying vector-resultant bending moment, which may only
29
resultant max mudline moment
Ratio of SRSS to vector− 1.2 360
1 180
0.8 0
0 60 120 180 240 300 360 0 60 120 180 240 300 360
Mean Wind Direction (degrees) Mean Wind Direction (degrees)
(a) (b)
Figure 2.10: (a) Ratio of the SRSS of the maximum mudline bending moments in two
orthogonal directions to the maximum vector-resultant mudline bending moment; and
(b) Variation of the mean direction of the vector-resultant mudline bending moment
with the mean wind direction.
be obtained from the campaign (i.e., time series) data. Figure 2.10a shows the ratio
of these two measures of the maximum mudline moment, and it is found that this
ratio is very close to unity, and that it rarely exceeds 1.1. Also, the mean direction of
the vector-resultant mudline bending moment is generally almost coincident with the
mean wind speed direction, as shown by Fig. 2.10b, and thus it can be assumed to
represent the fore-aft tower bending moment. Therefore, the SRSS of the maximum
mudline moments in the two orthogonal directions is a reasonable, albeit slightly con-
servative, approximation for the maximum of the vector-resultant mudline moment
and we use it as the load random variable in the present study and note that it closely
represents the maximum fore-aft mudline bending moment. Henceforth, we refer to
this load variable simply as the maximum mudline bending moment and denote it by
the symbol, M .
30
25 3
20
u (MN−m)
β (MN−m)
2
15
10
1
5
0 0
6 6
5 26 5 26
4 2224 4 2224
3 161820 3 161820
H 2 1214 H 2 1214
(m) 1 8 10 /s) (m) 1 8 10 /s)
s 0 4 6 V (m
s 0 4 6 V (m
(a) (b)
Figure 2.11: Estimates of the Gumbel distribution parameters, (a) u and (b) β, for
maximum mudline bending moment in the 90◦ -120◦ wind-direction sector.
25 3
20
u (MN−m)
β (MN−m)
2
15
10
1
5
0 0
6 6
5 26 5 26
4 2224 4 2224
3 161820 3 161820
H 2 1214 H 2 1214
(m) 1 8 10 /s) (m) 1 8 10 /s)
s 0 4 6 V (m
s 0 4 6 V (m
(a) (b)
Figure 2.12: Estimates of the Gumbel distribution parameters, (a) u and (b) β, for
maximum mudline bending moment in the 240◦ -270◦ wind-direction sector.
where u and β are parameters of the Gumbel distribution that represent the modal
value and a measure of dispersion, respectively. To establish these parameters, we bin
the data for each wind-direction sector according to a moving window representing
an interval on mean wind speed and on significant wave height. The moving window
31
Table 2.3: Distribution of datasets by wind-
direction sectors.
Number of (V, Hs ) bins
No. of that had n datasets
datasets, n 90◦ < θ < 120◦ 240◦ < θ < 270◦
0 521 648
0-20 128 76
20-40 58 70
40-60 44 6
60-80 19 0
80-100 12 0
100-120 7 0
120-140 5 0
140-160 2 0
160-180 3 0
180-200 1 0
0-200 800 800
for mean wind speed of 2 m/s is the same as was described earlier (and used in
Fig. 2.9). The window for significant wave height is of size 1 m; this window is moved
in increments of 0.25 m over the range from 0 m to 5 m. Thus, a total of 40 bins for
the mean wind speed random variable and 20 for the significant wave height random
variable are employed. For each bin, we estimate the parameters using least squares
regression with a transformed version of Eq. 2.4 carried out over the larger half of the
data (i.e., above the median fractile) so as to obtain a better representation of the
distribution tails since we are most interested in extrapolating from observed loads
to rarer unobserved extremes. As an illustration, Fig. 2.11 shows estimated values
of Gumbel parameters for one wind-direction sector (90◦ -120◦ ) for winds from the
sea and Fig. 2.12 shows estimates for one wind-direction sector (240◦ -270◦ ) for winds
from the shore. It is observed that, for both these wind-direction sectors, the modal
32
value, u, follows a smoothly increasing trend with increasing wind speed, which is
to be expected as turbine tower moments increase with wind speed. The Gumbel
dispersion parameter, β , on the other hand, fluctuates over a considerably wide
range of values; sometimes these estimates are large due to large dispersion in the
observed maximum mudline moment values in some bins, while at other times, they
are large because of extremely limited data in some bins. Table 2.3 summarizes the
distribution of available datasets over the different (V, Hs ) bins for the two wind-
direction sectors, (90◦ -120◦ ) and (240◦ -270◦ ). It can be clearly seen that for the
240◦ -270◦ wind-direction sector, no bins have more than sixty datasets, only 6 bins
have between forty and sixty datasets, and 146 bins have between zero and forty
datasets. Data distribution is somewhat better for the 90◦ -120◦ wind-direction sector
where several bins have more than 100 datasets. For both wind-direction sectors,
though, there are a large number of bins with no data at all. Such data distributions
as seen for the two sectors studied here, with slightly more data for wind-direction
sectors associated with winds blowing from the sea but with several (V, Hs ) bins with
sparse or no data both for winds from the shore and winds from the sea, is typical,
and serves to highlight the limited nature of the data.
33
Table 2.4: Statistics of the environment and load time series for the three
selected campaigns studied in Figs. 2.13, 2.14 and 2.15.
Wind Speed (m/s) Sea surface Extreme
Campaign Nacelle Met. Mast elevation MBM
Mean SD Mean SD Hs (m) Skew. (MN-m)
ShoreV18H0N1 19.1 4.0 15.5 3.6 0.64 -0.58 23.3
ShoreV18H0N9 18.6 2.4 13.9 2.2 0.64 0.17 14.9
SeaV12H3N41 13.2 1.8 12.7 1.7 3.47 -0.12 19.4
Note: SD: Standard deviation; Hs : Significant wave height; Skew.: Skew-
ness; MBM: Mudline bending moment.
j and (j + 1) m; the notation, Nk, represents the k th time series from this bin.
The prefixes “Shore” and “Sea” refer to winds from the shore and sea, respectively.
Statistics for the wind, waves, and mudline bending moment for these two time series
are summarized in Table 2.4 (one additional time series, referred to as ShoreV18H0N9,
is also discussed in the table for comparison purposes).
Table 2.4 shows that the ten-minute extreme mudline bending moments are
23.3 and 19.4 MN-m for the time series, ShoreV18H0N1 and SeaV12H3N41, re-
spectively. The largest load during the campaign was recorded for winds from the
shore. It is found that turbine loads are primarily influenced by wind speed and
only very slightly by wave height variation. It can also be seen from Table 2.4 that
the mean wind speeds recorded at the nacelle for time series, ShoreV18H0N1 and
SeaV12H3N41, are 19.1 and 13.2 m/s, respectively. Both these mean wind speeds
are above the rated wind speed of 12 m/s for this turbine. More interesting is the
difference in the standard deviations of the wind speed for these two time series. The
turbulence standard deviation for ShoreV18H0N1 is more than twice the value ob-
served for SeaV12H3N41. We point out here that wind speed measurements at the
34
nacelle are likely affected by the turbine wake. Only the mean wind speed from the na-
celle data was corrected to match the free-field mean wind speed [10]; the turbulence
standard deviation at the nacelle may not give an accurate estimate of the turbulence
character at hub height. Hence, we compare turbulence standard deviations from the
nacelle to that from another measurement of wind speed at a meteorological mast lo-
cated on the shore (Fig. 2.2a). These latter measurements were, however, at a height
of 40 m, which is lower than the hub height of 62 m, and also this met mast is about
1 km away from the turbine. Still, if either the nacelle or the met mast data are con-
sidered, the turbulence standard deviation in data segments for winds from the shore
is larger (by a factor of greater than two) than that in data segments for winds from
the sea. Since turbulence standard deviation (or equivalently, turbulence intensity,
which is the turbulence standard deviation divided by the mean wind speed) directly
influences turbine loads, we might expect that the largest loads might result during
conditions associated with winds from the shore and will likely occur for higher wind
speeds there.
Table 2.4 also shows the skewness of the sea surface elevation process. It may
be seen that all the campaigns have non-zero skewness values, which implies that
the irregular waves are asymmetric and hence nonlinear (and non-Gaussian); such
nonlinear nature of waves is characteristic of shallow water sites.
35
Wind Speed at Nacelle (m/s)
3 Wind Speed at Nacelle Wind
3 Speed at 40 m Met−mast
30 10 10
20
10 2 2
PSD, (m/s)2/Hz
PSD, (m/s)2/Hz
10 10
0
0 50 100 150 200 1 1
Wind Speed at 40 m Met−mast (m/s) 10 10
30 0 0
10 10
20
10 −1 −1
0 10 10
0 50 100 150 200 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Sea Surface Elevation, (m) Frequency (Hz) Frequency (Hz)
1 Sea Surface Elevation 2Mudline Bending Momemt
0.5 0.5 10
0
PSD, (MN−m)2/Hz
−0.5 0.4 1
−1 10
PSD, m2/Hz
0 50 100 150 200 0.3 0
Mudline Bending Momemt (MN−m) 10
0.2
20 −1
10
10 0.1
−2
0 0 10
0 50 100 150 200 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
time (sec) Frequency (Hz) Frequency (Hz)
Figure 2.13: Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the mudline bending
moment for the campaign time series, ShoreV18H0N1
PSD, (m/s)2/Hz
PSD, (m/s)2/Hz
10 10
0
0 50 100 150 200 1 1
Wind Speed at 40 m Met−mast (m/s) 10 10
30 0 0
10 10
20
10 −1 −1
0 10 10
0 50 100 150 200 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Sea Surface Elevation, (m) Frequency (Hz) Frequency (Hz)
1 Sea Surface Elevation 2Mudline Bending Momemt
0.5 0.5 10
0
PSD, (MN−m)2/Hz
−0.5 0.4 1
−1 10
PSD, m2/Hz
Figure 2.14: Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the mudline bending
moment for the campaign time series, ShoreV18H0N9
36
Figures 2.13, 2.14, and 2.15 show time series and power spectral density (PSD)
functions for the wind speed data (from both the nacelle and the met mast), the
sea surface elevation, and the mudline bending moment for the three selected time
series. Only a 200-second segment of the time series is shown where the largest
load was recorded in each case. As expected, the nacelle turbulence time series and
power spectra show more high-frequency content than is seen from the met mast
data; this is likely due to turbine wake effects. Comparing Figs. 2.13 and 2.14—both
for winds from the shore—it is seen that for these two time series with comparable
mean wind speeds, turbine loads are higher when the turbulence intensity is higher
(thus, ShoreV18H0N1 in Fig. 2.13 has higher loads than ShoreV18H0N9 in Fig. 2.14).
Both ShoreV18H0N1 and ShoreV18H0N9 have comparable significant wave heights
and both also have comparable sea surface elevation PSD peaks at the tower nat-
ural frequency (approximately, 0.5 Hz); clearly, since ShoreV18H0N1 has the larger
turbulence intensity, wind is the more important contributor to turbine loads than
waves. Still, resonant vibrations of the tower are at least somewhat contributed to
by waves and it is clear from studying Fig. 2.13, for example, that the largest load
in the mudline bending moment time series is recorded between 150 and 200 seconds
where the load process is oscillating at a frequency of around 0.5 Hz or a period of 2
seconds where waves provide some energy.
To summarize, for the two cases associated with winds from the shore, the
ShoreV18H0N1 case where the turbulence standard deviation is 3.6 m/s (at the met
mast) experiences the larger load of 23.3 MN-m compared to 14.9 MN-m in the
V18H0N9 case where the turbulence standard deviation is 2.2 m/s.
Figure 2.15 shows similar time series and PSD estimates for a data segment
associated with winds from the sea, SeaV12H3N41, as were studied in Figs. 2.13
37
Wind Speed at Nacelle (m/s)
30
2 Wind Speed at Nacelle Wind
2 Speed at 40 m Met−mast
20 10 10
10 1 1
PSD, (m/s)2/Hz
PSD, (m/s)2/Hz
0 10 10
0 50 100 150 200
Wind Speed at 40 m Met−mast (m/s) 0 0
10 10
30
20 −1 −1
10 10
10
−2 −2
0 10 10
0 50 100 150 200 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Sea Surface Elevation, (m)
Frequency (Hz) Frequency (Hz)
4
2 Sea Surface Elevation 2Mudline Bending Momemt
15 10
0
−2
PSD, (MN−m)2/Hz
1
−4 10
PSD, m2/Hz
0 50 100 150 200 10
Mudline Bending Momemt (MN−m) 0
10
20 5 −1
10
10
0 −2
0 10
0 50 100 150 200 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
time (sec) Frequency (Hz) Frequency (Hz)
Figure 2.15: Time series and power spectral density (PSD) functions of the wind
speed at the nacelle and met mast, the sea surface elevation, and the mudline bending
moment for the campaign time series, SeaV12H3N41
and 2.14 for the cases with winds from the shore. The largest recorded load for
winds from the sea occurred for the selected data segment. As was discussed earlier,
winds from the sea are less turbulent than winds from the shore; here, for example,
the turbulence standard deviation is only 1.7 m/s at the met mast compared to the
values of 3.6 m/s and 2.2 m/s, respectively, for the cases for winds from the shore,
ShoreV18H0N1 and ShoreV18H0N9. As a result, despite the fact that the significant
wave height for SeaV12H3N41 is considerably larger than for ShoreV18H0N1, the
extreme turbine loads are smaller for this case (19.4 MN-m) compared to 23.3 MN-m
for ShoreV18H0N1. Again, we reiterate that turbulence has primary influence on
turbine loads and waves have only secondary influence.
While only three ten-minute time series were selected here to illustrate differ-
ences between winds from the sea and the shore and their influence on turbine loads,
the conclusions reached from studying these—namely that winds from the shore cause
larger loads due to their higher turbulence intensities and that wind speed is more
38
important than wave height—are supported when one studies other data segments
in the campaign as well. Consequently, we might expect that the long-term loads,
derived from statistical extrapolation, would also be larger for winds from the shore.
−4
10
1 year
−6 20 year
10
50 year
−8
10
10 20 30 40 50
Mudline Bending Moment (MN−m)
Figure 2.16: Comparison of probability of load exceedance curves for winds from the
sea, winds from the shore, and winds from all directions.
39
loads based on the “all-direction” winds are almost the same as those based on winds
from the shore. This is because at longer return periods, the exceedance probability
(of specified load levels) for winds from the shore is about two orders of magnitude
higher than that for winds from the sea. Thus, the exceedance probability for the
all-direction winds that represents the sum of these probabilities for winds from the
sea and the shore, is almost same as that for the winds from the shore alone.
While the analysis of this particular data set has shown that winds from the
shore govern the long-term loads (mudline bending moment), we are also interested in
determining which wind-direction sector accounts for these derived loads. Figure 2.17
shows long-term loads for return periods of 1 year and 20 years. It can be observed
that the relative importance of the different wind-direction sectors is generally the
same for both return periods. Focusing only on the 20-year return period, we find that
this load, if the 240◦ -270◦ wind-direction sector alone were representative, would be
about 36.5 MN-m, while that for the 180◦ -210◦ wind-direction sector would be about
35 MN-m. The actual 20-year load of 38 MN-m is slightly higher when considering
all winds from the shore (as well as considering all-direction winds, from the sea as
well as from the shore). Also, the 20-year load considering winds from the sea alone
is only about 30 MN-m, which is roughly 20% smaller than the 20-year load for winds
from the shore. From these results, we can conclude that the 240◦ -270◦ wind-direction
sector associated with winds from the shore governs the overall 20-year loads for this
turbine.
Since the parameters for the short-term load distribution (according to Eq. 2.4)
are statistically estimated from limited loads data, there is uncertainty associated
40
Figure 2.17: Long-term loads (mudline bending moments) for 1-year and 20-year
return periods for different wind-direction sectors. Also shown are 90%-confidence
intervals along with mean values for the ‘with-uncertainty’ case based on bootstrap
methods applied to the limited data in each sector.
with these estimates. This parameter uncertainty must clearly result in variability
in the predicted long-term loads as well. In an effort to quantify this variability, we
employ non-parametric bootstrap techniques [20] that rely on random resampling of
the data. We estimate parameters for the short-term load distribution from each of
Nb separate resamplings and, using this set of Nb (taken to be 100 here) estimates,
it is then possible to obtain statistics such as the standard error in each parameter
estimate that can help to quantify the variability in that parameter. Variability in
short-term distribution parameters propagates to variability in the long-term load
exceedance distribution. This can be quantified by computing load exceedance prob-
ability curves for each of the resamplings, from which we can then compute, say,
the 90-percent confidence interval on a long-term load associated with any specified
return period. Note that this estimated uncertainty in long-term loads as derived
41
from such a bootstrap analysis does not account for a more fundamental source of
uncertainty that might arise due to a severely limited data set, such as the use of only
2,300 data sets recorded over a total of about sixteen months, especially since multiple
harsh-weather seasons are not covered then. Besides, other simplifying assumptions
such as that of the use of the SRSS rule for combining orthogonal maximum bend-
ing moment components, for example, introduce additional uncertainty which is not
explicitly accounted for.
42
samplings. We use least-squares regression over the larger half of each loads data
resampling, which relies then, intentionally, on extreme values rather than on the
body of the overall distribution. When loads data are resampled, parameters based
on extreme values exhibit relatively large variability from one resampling to the next.
This problem becomes especially acute when the number of available data is small,
as is the case the data in some bins. To investigate this issue further, we estimate
parameters of the Gumbel distribution using three additional methods (to the first),
and recompute long-term loads and associated uncertainty bounds. The four parame-
ter estimation methods include: (a) regression over the above-median loads data (this
method is what we have employed in the preceding discussions); (b) regression over
all the loads data; (c) a collocation method that only uses 50th and 90th load per-
centiles; and (d) the method of moments (that is based on all the data). While details
related to each parameter estimation method may be obvious from the descriptions
of these methods, it is worth noting that the collocation method based on 50th and
90th percentile load levels effectively involves estimating parameters by drawing a line
through the 0.5 and 0.9 load fractiles on Gumbel probability paper. This yields the
Gumbel parameters, u and β, as follows:
m90 − m50
β= ; u = m90 + βp90 (2.5)
p90 − p50
³ ³ n ´´
pn = log − log ; n = 50, 90
100
where m50 and m90 are the median and 90-percentile levels of the mudline bending
moment, respectively, and pn is a function of the n-percentile load level. Note that
estimation methods (a) and (c) rely more heavily on the tails or extremes of the loads
while methods (b) and (d) make use of the full body of the distribution on loads.
Figure 2.18 shows a comparison of 90-percent confidence intervals (or uncertainty
43
bounds) on the 20-year loads based on these four methods. As was mentioned pre-
viously, uncertainty bounds for almost all the wind-direction sectors are significantly
larger when based on methods that emphasize on load extremes (methods (a) and (c)
here), compared to those based on methods that utilize all the data in the body and
the tails of the distribution. This is consistent with our observation that coefficient of
variation estimates (not shown here) on distribution parameters obtained from meth-
ods focused on load extremes are about 1.5-2 times larger than those obtained from
methods that rely on all the loads data. Note that mean estimates of the long-term
loads and the relative importance of the different wind-direction sectors do not vary
very much when the different estimation methods are used.
44
Figure 2.18: Comparison of 90-percent confidence intervals and mean value estimates
of the 20-year load based on different methods of estimation of parameters of the short-
term load distribution: (a) regression over the above-median loads data; (b) regression
over all the loads data; (c) collocation using the 50th and 90th load percentiles; and
(d) the method of moments.
45
2.5.2 Importance of Winds from the Shore
We now discuss possible physical explanations for why winds from the shore
govern long-term loads for this turbine. While the turbine in question is subjected
to both wind and wave loads, previous studies [10, 33] of the same Blyth data have
suggested that the contribution from waves to the overall tower loads is relatively
small. While these cited studies recognize that large wave heights are likely at the
Blyth site, they do not account for the possibility of breaking waves that might
explain large turbine loads. Nevertheless, we assume based on these studies that
wind loads at this site are more important as far as the mudline bending moment is
concerned. We thus focus our attention on the wind loading. The data show that
the mean turbulence intensity (measured at the 40-meter elevation station of the
meteorological mast, where measurements are not affected by the turbine wake) for
winds from the shore is about 0.21, while that for winds from the sea is about 0.09.
We also showed, while discussing several representative time series in Table 2.4, that
large loads often occur when the associated winds are more turbulent, as is the case
for winds from the shore. Therefore, the larger turbulence intensity accompanying
winds from the shore may be partly responsible for the larger long-term turbine loads.
46
are governed mainly by characteristics and topography of the land, resulting in high
turbulence. While the effect of turbulence on long-term loads for wind turbines sited
on land has been rather well studied (see, for example, [68]), we find here that long-
term loads for an offshore wind turbine especially when sited close to the shoreline
may also be significantly influenced by turbulence.
47
Figure 2.19: Comparison of 20-year loads when the environmental random variable
vector for each wind-direction sector, X, consists of (a) mean wind speed (V ) and
significant wave height (Hs ); and (b) mean wind speed (V ) and turbulence standard
deviation (σV ).
10%, which is very small. This suggests that turbulence standard deviation does not
have as significant an influence on turbine long-term loads as does mean wind speed.
For the mudline bending moment, then, we can say that the important environmental
random variable is primarily the mean wind speed, V , and even though there is some
influence of turbulence standard deviation, σV , this effect and that of waves are only
of secondary importance.
Our objective in this chapter was to derive long-term loads for a 2MW offshore
wind turbine sited in 9 meters of water at the Blyth site in the United Kingdom. Field
data were available on wind speed, wave height, and mudline bending moment in two
formats, as summary data (ten-minute statistics) and as campaign data (ten-minute
48
time series). The data were broadly classified into winds from the sea and winds
from the shore. Furthermore, the summary data were also divided into 30◦ wind-
direction sectors. Short-term load distributions conditional on environmental random
variables were established using ten-minute global maxima that were available from
the summary data. Integration of these distributions with the relative likelihood
of different wind speed and wave height combinations allowed derivation of turbine
long-term loads.
• Winds from the shore govern the long-term loads for the instrumented wind
turbine at Blyth. A detailed analysis of the available time series shows that
winds from the shore result in larger loads, compared to winds from the sea, due
to: (1) the larger turbulence intensity associated with winds from the shore; and
(2) the relatively greater amounts of energy in the waves at the tower natural
frequency (as seen, for example, by studying load time series associated with
winds from the shore compared to ones associated with winds from the sea).
49
In the design of wind turbines, statistical extrapolation is needed to predict
long-term loads according to design guidelines [16, 42, 43]. The loads data needed for
such extrapolation are obtained using stochastic turbine response simulations (and
not from field measurements). In light of conclusions reached from analysis of the
Blyth field data, the following questions need careful study when long-term turbine
loads are to be predicted using simulations.
• Most offshore wind farms are sited in fairly shallow waters. Waves at such sites
are nonlinear and need also to be represented as irregular (stochastic) processes.
Therefore, appropriate irregular nonlinear wave models should be employed to
compute hydrodynamic loads in turbine response simulations.
Both of these issues are addressed in the following chapters of this dissertation.
50
Chapter 3
3.1 Introduction
51
which type of load extremes—e.g., global maxima or block maxima (described in
Section 1.4.1)—should be extracted from each simulated ten-minute time series; an
additional question with use of block maxima relates to the choice of an appropriate
block size. Furthermore, it is of interest to know how these extreme models are
related and whether long-term loads predicted by different models for extremes are
comparable or not.
In this chapter, we use the LE3 data set to address some concerns that have
been raised with regard to experiences practitioners have had with attempts to address
DLC 1.1 in the IEC standard. In particular, we address efficiency in simulations when
we discuss which wind speeds are design drivers for various loads on a wind turbine; we
address convergence criteria that lead to approaches to quantitatively identify when
an adequate number of simulations have been run that can yield stable short-term
empirical load distributions; and we address the issue of independence in block sizes
52
using a statistical test for independence and also discuss differences in load predictions
based on the use of global versus block maxima. Throughout, insights are provided to
guide the effort involved in carrying out statistical loads extrapolation as required for
DLC 1.1. All of the findings that are reported here are based exclusively on statistical
studies on data from simulations with the baseline 5MW turbine model.
Design Load Case (DLC) 1.1 in the wind turbine standard IEC 61400-1 [42]
deals with loads for ultimate limit states, to be estimated using extrapolation, that
a wind turbine might experience during normal operation—when wind speeds range
between cut-in and cut-out, the minimum and maximum wind speeds, respectively,
at which a wind turbine operates. DLC 1.1 also specifies that inflow conditions used
in the simulations should represent those associated with near-neutral atmospheric
conditions. A normal turbulence model (NTM) is prescribed in the IEC standard
that must be used to represent the inflow turbulence fields pertinent to DLC 1.1.
In this load case, the hub-height wind speed, V , averaged over ten minutes may be
treated as a single random variable representing the environment; thus the vector
of environmental random variables, X, described in Eq. 1.1, consists of only one
variable, V . DLC 1.1 requires that aeroelastic simulations be conducted over the
entire normal operating wind speed range. It is convenient, as the IEC standard
permits, to carry out simulations over discrete wind speed intervals or bins; typically,
bins of 2 m/s size for V are employed. Next, we discuss models for simulation of the
wind and of the turbine response.
53
3.2.2 Simulation Model for Inflow Wind
We briefly discuss the model employed for simulation, according to DLC 1.1,
of the three-dimensional stochastic wind velocity field on the turbine’s rotor plane,
which is normal to the longitudinal wind direction. The stochastic wind velocity field
is simulated on a grid of points on the rotor plane and, at each of these points, the
wind velocity has components in three orthogonal directions — longitudinal, lateral,
and vertical. Such a wind field is a multivariate stochastic process, which can be
simulated using the power spectrum, Sk (f ), of each component of wind velocity at a
specified location, where the subscript, k, denotes the direction (1 = longitudinal, 2
= lateral, and 3 = vertical) and f represents frequency in Hertz. The wind velocity at
any two points in space is correlated and, therefore, we also need frequency-dependent
coherence functions for the turbulence components and for different spatial separa-
tions between points on the rotor plane. Such a multivariate stochastic process can
be simulated using inverse Fast Fourier Transform techniques; algorithms for such
simulation procedures are readily available in the literature [94]. A detailed proce-
dure for simulation of the stochastic wind field for wind turbines is also described by
Saranyasoontorn [87], and we do not repeat it here. We limit our presentation here
to introducing specific models, as suggested by the IEC standard, that are commonly
used for simulation of the wind field.
We use a Kaimal spectrum [49], per the IEC standard, to describe the one-sided
power spectral density functions for each of the three components of wind velocity.
The spectrum in normalized form is given as follows.
f Sk (f ) 4f Lk /V
= (3.1)
σk2 (1 + 6f Lk V )5/3
where V is the mean wind speed at hub height. The variables, σk and Lk , are the
standard deviation and the integral scale parameter, respectively, of the k th wind
54
velocity component. Also, σk2 represents the area under the power spectrum, Sk (f );
for the Kaimal model, σ22 = 0.64σ12 , while σ32 = 0.25σ12 , while σ1 is defined according
to the wind turbine class and reference turbulence intensity [42]. The integral scale
parameter, Lk , is related to the longitudinal turbulence scale parameter Λ1 , such that
L1 = 8.1Λ1 , L2 = 2.7Λ1 , L1 = 0.7Λ1 . The longitudinal turbulence scale parameter,
Λ1 , at height z, is given as follows.
½
0.7z z ≤ 60 m
Λ1 = (3.2)
42 m z > 60 m
The turbulence scale parameter represents the size of turbulent eddies, and it depends
on the height above ground level and also on the surface roughness if the height at
which the wind velocity is being simulated is close to the ground [9].
where Lc , the coherence scale parameter, is the same as the integral scale parameter,
L, for the longitudinal component. Also, r refers to the separation between two points
projected on a plane normal to the average wind direction.
The mean wind speed in the longitudinal direction varies with height, z; such
a variation of the mean wind speed with height, V (z), can be described by the fol-
lowing power law, which is referred to as the ‘normal wind profile model’ in the IEC
standard [42]. Thus, we have
³ z ´α
V (z) = V (3.4)
H
where V is the mean wind speed at the hub height, H, and the power law exponent,
α, is assumed to be 0.2.
55
The turbulence standard deviation of the longitudinal wind velocity compo-
nent, σ1 , is a random variable; a lognormal model can be assumed for the probability
distribution of σ1 conditional on the mean wind speed V [42]. However, the IEC stan-
dard permits the use of a representative turbulence standard deviation corresponding
to the 90-percentile value of σ given V . This model is referred as the ‘Normal Tur-
bulence Model’ (NTM) and the representative σ1 value in m/s is given as follows:
where Iref is the expected value of the turbulence intensity (turbulence standard
deviation divided by the mean wind speed) at a reference mean wind speed of 15
m/s; also b = 5.6 m/s and V is the hub-height mean wind speed in m/s . The IEC
standard defines three categories of turbulence, namely A, B, and C. The values of
Iref are 16%, 14% and 12% for turbulence categories, A, B, and C, respectively.
The spectral model for the inflow wind velocity field on a wind turbine rotor
and the procedure for its simulation in the time domain are built into several computer
programs. One such program, TurbSim [45, 51], developed at the National Renewable
Energy Laboratory, was used to simulate the stochastic wind velocity field for the
LE3 loads data sets. TurbSim supports all the features of an earlier program for
simulation of wind, SNLWIND-3D [50], which itself was based on the original program
SNLWIND [104].
56
wind turbines—both onshore and offshore—in the time domain. FAST is a medium-
complexity code that models a wind turbine as a combination of rigid and flexible
bodies wherein flexible beam elements (representing the blades and the tower) are
formulated in terms of generalized coordinates. FAST models various control sys-
tems including blade and pitch control, yaw control, and high-speed shaft brake
control. Another program that may be used for wind turbine response simulation
is ADAMSr , developed by MSC Corporation, which can perform detailed finite el-
ement analysis of turbine elements such as blades. Other wind turbine response
simulators include BLADED [7] from Garrad Hassan and Partners, FLEX5 [77] from
the Technical University of Denmark, HAWC2 [54] from Riso National Laboratories,
and PHATAS [56] from the Energy Research Centre of the Netherlands (ECN). The
LE3 group used FAST for wind turbine response simulations; later in this disserta-
tion, we also use FAST for our turbine response simulations. FAST is open-source
software and its source code as well as an executable file are available from the NREL
web site (http://wind.nrel.gov/designcodes/simulators/fast).
The structural model for a wind turbine in FAST consists of five flexible bodies
(the tower, three blades, and a drive shaft) and up to nine rigid bodies (earth, support
platform, base plate, nacelle, armature, gears, hub, tail and structure furling with the
rotor). These flexible and rigid bodies together define 24 degrees of freedom (DOFs),
based on a combined modal and multibody dynamics formulation. The translational
(surge, sway, and heave) and rotational (roll, pitch, and yaw) motions of the support
platform, relative to the inertial reference frame (earth), describe the first six DOFs.
The first two fore-aft (longitudinal) and the first two side-to-side (lateral) bending
modes of the tower describe the next four DOFs. The yawing motion of the nacelle,
the generator azimuth angle, and the drivetrain compliance (between generator and
57
rotor) provide the next three DOFs. Each of the three blades is modeled by its first
two flapwise bending modes and its first edgewise bending mode, providing a total
of nine DOFs for the three blades. The rotor-furl and tail-furl provide the last two
DOFs, although these DOFs are typically not modeled for modern wind turbines.
The FAST program allows a user to choose the necessary combination of DOFs to
appropriately model any horizontal-axis wind turbine.
An integral part of modern wind turbines are the control systems, which are
mainly designed to regulate the rotor speed in order to optimize power production
from the generator. One of the most common control systems is the blade pitch
control system, in which, over some wind speed ranges, the pitch angle of a blade
is increased in response to an increase in the instantaneous wind speed, so as to
reduce aerodynamic loads on the blades and consequently reduce the rotor speed.
Other forms of control strategies that can be modeled in FAST include control of the
generator torque, application of brakes on the high speed shaft, deployment of brakes
on blade tips, nacelle yaw. Each wind turbine can have a different control system that
can be algorithmically implemented in FAST by including appropriate subroutines.
58
of this theory can be found in Saranyasoontorn [87]. Once these forces are available,
FAST uses numerical techniques to solve the dynamic equations of motion for the
wind turbine system in the time domain to yield the turbine response.
The Loads Extrapolation Evaluation Exercise (LE3 ) working group was consti-
tuted in order to address ongoing issues related to the IEC standard and particularly
DLC 1.1 that deals with statistical loads extrapolation from limited simulation. The
LE3 working group approved a proposal to develop a database of simulated load time
histories and summary statistics from a 5MW turbine model developed by the Na-
tional Renewable Energy Laboratory (NREL). This model is based on an onshore
version of NREL’s baseline turbine model developed to represent a utility-scale 5MW
offshore wind turbine [48]; this onshore model has identical properties to the offshore
turbine model (extensively used in later chapter of this dissertation) above the mean
sea level. The turbine has a hub height, above the ground, of 90 meters and a rotor
diameter of 126 meters. The machine is a variable-speed, collective pitch-controlled
turbine with a rated wind speed of 11.5 m/s. The maximum rotor speed is 12.1 rpm.
Moriarty [67] provides a detailed description of the turbine model and the various
inflow conditions covered by the simulations. Inflow turbulence was simulated using
TurbSim; along with the wind model discussed in Section 3.2.2. The program, FAST,
was used to carry out aeroelastic simulations for hub-height mean wind speeds, V ,
varying between cut-in and cut-out mean wind speeds. For the IEC Class I-B site
assumed, the ten-minute hub-height mean wind speed follows a Rayleigh distribution
with average equal to 10 m/s.
The data set used in this study, which was generated by the LE3 working group,
59
consists of 14,400 ten-minute simulations or 1,200 simulations for each of twelve wind
speeds ranging between 3 and 27 m/s. When running TurbSim to generate inflow
turbulence time histories, the target wind speeds were set at discrete values of 3 m/s,
5 m/s, etc. up to 25 m/s and were assumed to represent 2 m/s bins centered at
the target values. Realized ten-minute average wind speeds varied slightly from the
target values. The simulations were run at NREL using a distributed computing grid
consisting of sixty desktop computers, such that a single ten-minute simulation was
completed, on average, in about 10 seconds. Note that if only one desktop computer,
of today’s typical configuration, is used, it takes about 8 minutes to complete one
ten-minute simulation.
In the present study, four representative and contrasting loads were ana-
lyzed from the LE3 data set. These loads include the out-of-plane bending moment
(OoPBM) at a blade root, the out-of-plane blade tip deflection (OoPTD), the fore-aft
tower bending moment (FATBM) at the base, and the in-plane blade bending mo-
ment (IPBM) at a blade root; we will use the acronyms to describe these four load
types throughout this chapter.
To compare the relative importance of different wind speed bins on load ex-
tremes, it is of interest to study load extreme statistics as a function of wind speed.
With a little effort (i.e., limited simulations), it is often possible to identify which
wind speeds cause the largest turbine loads on average and, equally important, which
ones show the greatest load variability. The largest loads are associated with the low-
est probability of exceedance and, as such, are closest to the rare probability levels to
which extrapolation is needed; the wind speeds where these largest loads occur are
60
20 100
OoPBM (MN−m) 80
FATBM (MN−m)
15
60
10
40
5 Mean Mean
Max 20 Max
Min Min
0 0
2 4 6 8 10 12 14 16 18 20 22 24 26 2 4 6 8 10 12 14 16 18 20 22 24 26
Mean Wind Speed (m/s) Mean Wind Speed (m/s)
10 10
8
8
IPBM (MN−m)
OoPTD (m)
6
6
4
4 Mean Mean
Max 2 Max
Min Min
2 0
2 4 6 8 10 12 14 16 18 20 22 24 26 2 4 6 8 10 12 14 16 18 20 22 24 26
Mean Wind Speed (m/s) Mean Wind Speed (m/s)
Figure 3.1: Distribution of short-term load maxima as a function of wind speed for
four loads, FATBM, OoPBM, IPBM and OoPTD, based on 200 simulations per wind
speed bin.
From Fig. 3.1, it is possible to identify those wind speed bins associated with
the largest loads and greatest variability in short-term load extremes. This is useful
since summary plots such as these make it possible to focus efforts in the most im-
portant bins. For OoPBM and FATBM, important controlling wind speed bins are in
61
the 14-22 m/s range with perhaps the dominant winds being closest to the 14-16 m/s
bin. The OoPTD is controlled by somewhat lower wind speeds; largest loads occur
between 10 and 16 m/s. Note that the behavior of a pitch-controlled wind turbine,
such as the one studied here, changes when the inflow wind speed exceeds the rated
wind speed, i.e., the wind speed at which the blade-pitch control action start to act;
the rated wind speed for our 5MW turbine is 11.5 m/s, and thus it can be inferred
that wind speeds near and above the rated wind speed are important for loads such as
OoPBM, FATBM and OoPTD. The IPBM, on the other hand, is clearly dominated
by wind speeds close to the cut-out wind speed of 25 m/s. While these controlling
wind speed bins were identified by carrying out 200 simulations per bin, it is possible
to identify important bins with considerably less simulation effort. From Fig. 3.1, it
is clear that wind speeds below 10 m/s do not contribute large loads to any of the
four load types discussed; moreover, they also do not exhibit large variability in load
extremes; therefore, these low wind speed bins are not important in extrapolation
and a minimal effort in terms of simulation may be justified for these bins when de-
riving short-term distributions based on Eq. 1.1. In carrying out simulations with a
view towards extrapolation, it is worthwhile to understand turbine load extremes as
a function of wind speed in this manner to avoid excessive computational effort.
An alternative to the use of only a single maximum (i.e., the global maximum)
from each simulated ten-minute time series is to extract several extremes from each
time series in a systematic manner. While this can be done by methods such as the
peak-over-threshold procedure [83], there are simpler methods as well. For instance,
one could split or partition the time series into individual non-overlapping blocks
62
14
12
OoPBM (MN−m)
10
6
Block Maxima Global Maxima
4
1 2 3 4 5 6 7 8 9 10
time (minutes)
Figure 3.2: Example OoPBM load time series showing global and block maxima.
of constant duration [15]. From each of these blocks, then, a single largest value is
extracted; together these extracted extremes constitute a set of block maxima. Fig-
ure 3.2 shows an example load time series for out-of-plane bending moment (OoPBM)
with one-minute block maxima and the single ten-minute global maximum.
Figure 3.2 indicates that from a single load time series of ten-minute duration,
more extremes data may be extracted if block maxima are employed as opposed to
global maxima in extrapolation. As long as the block maxima can be shown to be
mutually independent, the short-term distribution, FL|V =v (l), of global maxima, L,
can be related to the short-term distribution, FLBlock |V =v (l), of block maxima, Lblock
as follows:
£ ¤n
FL|V =v (l) = FLBlock |V =v (l) (3.6)
63
0 OoPBM; 16 < V < 18 m/s; 20 Simulations; Block size = 20 sec
10
Block Maxima
Global Maxima
−1
10
Exceedance Probability −2
10
p80 − Block Max
−3
10
−4
10
5 10 15 20
OoPBM (MN−m)
Figure 3.3: Comparison of block and global maximum probability levels associated
with a given load quantile for OoPBM loads.
64
represent global maxima and the dots, block maxima. To highlight the point that
the global maxima are also block maxima, the specific global maxima extracted are
shown twice to indicate where they appear in the block maxima distribution. Some
extracted block maxima are higher than a few global maxima and arguably better
defined tail trends are seen in the block maxima distribution. However, as can be
seen, the 80th percentile ten-minute maximum load is at an exceedance probability
level of 1 − 0.978 or 0.022 if the block maxima distribution is used.
As was stated before, Eqs. 3.6 and 3.7 are valid as long as block maxima
selected from each time series are independent of each other. Intuitively, it may be
expected that smaller block sizes will lead to greater dependence among the block
maxima. Statistical tests for independence represent the only objective means of
assessing the extent of independence or lack thereof in a sample of block maxima
from load simulations.
Several tests to evaluate independence between two random variables are avail-
able in the literature. We use a test proposed by Blum et al. [5], which is a generalized
version of the test proposed by Hoeffding [34]. Details related to this test along with
examples may be found in Hollander and Wolfe [36]; that reference also provides a
correction for a typographical error in an equation in Blum et al [5]. Blum’s test has
been used by Skaug and Tjøstheim [95] to test for independence in time series data,
for which it was not originally developed.
65
sis, H0 , is that the two variables, X and Y , are independent. Thus, we have in terms
of cumulative distribution functions:
Blum’s test makes use of a test statistic, B, that must be checked against a
critical value, Bcr , at any specified significance level. This test statistic is computed
as follows:
N
X (N1 N4 − N2 N3 )2
1
B = π4N (3.9)
2 j=1
N5
where N is the sample size for both X and Y . The quantities, N1 to N4 , are computed
for all values of j from 1 to N or effectively for all choices of (X, Y ) = (xj , yj ) such
that
If the B value as computed from Eq. 3.9 is higher than Bcr , then the null hypothesis
is rejected and the two variables, X and Y , are not independent at the specified
significance level.
Useful illustrative examples of the use of this test involve examining this B
statistic for paired data that are known to be either strongly dependent or indepen-
dent. Note that for a bivariate Gaussian distribution, correlation implies dependence
(the distribution is completely defined by a correlation coefficient and the first two
66
N1 = 22; N4 = 22; N2 = 4; N3 = 2 N1 = 11; N4 = 21; N2 = 8; N3 = 10
3 2
1.5
2
N3 N4 1 N3 N4
1
0.5
(xj,yj)
0 0 (xj,yj)
y
y
−0.5
−1
N1 N2 −1 N1 N2
−2
−1.5
−3 −2
−3 −2 −1 0 1 2 3 −2 −1 0 1 2
x x
(a) (b)
Figure 3.4: Scatter plot of simulated samples of two bivariate Gaussian random vari-
ables with correlation coefficients of (a) 0.9 and (b) 0.05. Also indicated are the values
of N1 , N2 , N3 , and N4 as computed when carrying out Blum’s test for independence.
marginal moments of each variable). If Blum’s test is carried out for two jointly
distributed Gaussian random variables that are strongly correlated, the B statistic is
likely to be large; the opposite is true if the correlation is weak.
67
OoPBM FATBM
20 20
10<V<12 10<V<12
16<V<18 16<V<18
15 22<V<24 15 22<V<24
B
B
10 10
5 1% Level 5 1% Level
0 0
20 40 60 80 100 120 20 40 60 80 100 120
Block size (sec) Block size (sec)
IPBM OoPTD
20 20
10<V<12 10<V<12
16<V<18 16<V<18
15 22<V<24 15 22<V<24
B
10 B 10
5 1% Level 5 1% Level
0 0
20 40 60 80 100 120 20 40 60 80 100 120
Block size (sec) Block size (sec)
Figure 3.5: Variation of Blum’s test B statistic for four loads as a function of block
size (computed from 200 ten-minute time series for each load type and in three wind
speed bins).
computed B value is relatively smaller than in the first case. For this case, B is equal
to 2.80, which is smaller than the critical value, Bcr , of 4.23 at a 1% significance level.
Hence, the null hypothesis of independence is not rejected. These two examples serve
to illustrate the use of Blum’s test for independence between two random variables
in general.
The independence of block maxima of turbine loads for different block sizes
may be studied in a similar manner to that used for the preceding illustrative ex-
68
Table 3.1: Suggested block sizes (in seconds) for independent block maxima based
on mean (µB ) and mean plus one standard deviation (µB + σB ) values from 200
simulations and tested at the 1% significance level.
Wind speed OoPBM FATBM IPBM OoPTD
(m/s) µB (µB +σB ) µB (µB +σB ) µB (µB +σB ) µB (µB +σB )
2<V <4 50 70 30 50 30 60 50 70
4<V <6 40 60 25 40 40 60 40 60
6<V <8 40 60 30 50 30 50 40 60
8 < V < 10 40 50 40 50 25 40 40 50
10 < V < 12 15 20 20 25 20 30 15 20
12 < V < 14 20 30 20 30 15 20 25 40
14 < V < 16 20 30 20 30 15 20 20 30
16 < V < 18 15 25 15 25 15 20 20 25
18 < V < 20 15 20 6 15 10 15 15 25
20 < V < 22 7 20 5 6 9 15 15 20
22 < V < 24 7 15 4 5 8 15 15 20
24 < V < 26 6 15 4 5 7 15 15 20
ample. Blum’s test statistic, the B value, for block maxima may be computed by
forming lag-one vectors, X and Y , from all the block maxima in each ten-minute time
series. The B value may be computed for these lag-one extremes to test if they are
independent. It is expected that these extremes will become more independent as
the block size is increased. At a certain optimum block size, computed B values will
fall below the critical value, Bcr . While it is possible to study the B values for each
simulation corresponding to a given wind speed, it is more instructive to study these
B values (and, thus, independence) statistically as a function of block size by consid-
ering multiple simulations at a time in each wind speed bin; this makes it possible to
consider scatter or uncertainty in the B values over different simulations. To this end,
the mean, µB , and standard deviation, σB , of the B values from 200 simulations are
computed for four different load measures, OoPBM, FATBM, OoPTD, and IPBM.
69
Note that even with a small number of simulations, on the order of 15 to 20, statistics
of the B values are quite stable and 200 simulations are not really needed. The mean
B values with error bars representing one standard deviation are shown in Fig. 3.5
for the four load types and for three different wind speed bins, 10-12 m/s, 16-18
m/s, and 22-24 m/s. As expected, mean values of B decrease monotonically with
increasing block size. Even if the more stringent (µB + σB ) level is checked against
the critical value, Bcr , at the 1% significance level, independence of block maxima is
virtually assured for block sizes longer than 30 seconds for all four load types and in
all three wind speed bins. Summarized in Table 3.1 are the appropriate block sizes
based on criteria where either µB or (µB + σB ) values are compared to Bcr at the 1%
significance level. Clearly, loads in some wind speed bins exhibit greater dependence
than others (e.g., low wind speed bins) but it appears that, at least with the LE3
loads data set, one could safely choose block sizes of around 40-60 seconds, extracting
between 10 and 15 extremes (block maxima) from each ten-minute time series and
using them to establish short-term load distributions.
70
3.5.3 Discussion on Independence
It was stated earlier that independence of block maxima is required for the
mathematical relationships expressed in Eq. 3.6 to hold. However, it is also important
and perhaps more useful to discuss what effect, if any, the assumption of independence
(whether or not it is justified) has on the rare load fractiles that are of interest for
extrapolation. One can study, for example, the effect of assuming independence, or
of assuming a block size that assures independence, on the tails of short-term load
distributions. Table 3.3 shows estimates of the 84th percentile global maximum load
for four load types obtained using block maxima distributions for different block sizes,
while adjusting the fractile level according to Eq. 3.7.
Table 3.3 reveals that the selection of block size has almost no effect on short-
term loads at fairly rare fractile levels. This suggests that while smaller block sizes do
not guarantee independence, they do not predict grossly inaccurate loads, at least at
the 84th non-exceedance probability level for global maxima. This suggests that one
could use very small block sizes without concern for independence but no new informa-
71
Table 3.3: Estimates of the 84th percentile global maximum load for four different
load types as obtained from block maxima distributions using different block sizes.
Block OoPBM FATBM IPBM OoPTD
size (sec) 16 < V < 18 16 < V < 18 22 < V < 24 16 < V < 18
5 13.59 81.50 7.61 7.53
10 13.59 81.50 7.61 7.53
15 13.59 81.50 7.61 7.53
20 13.59 81.50 7.61 7.53
30 13.59 81.51 7.61 7.53
60 13.60 81.51 7.61 7.53
tion is gained by the significantly larger sample of extremes that would result with this
small block size. Stated differently, the findings summarized in Table 3.3, based on
this LE3 data set, suggest somewhat surprisingly that enforcing independence among
block maxima will likely have no significant impact on short-term distributions and
hence on extrapolation.
0 16 < V < 18 m/s; Block size = 20 sec 0 16 < V < 18 m/s; Block size = 5 sec
10 10
Block Maxima Block Maxima
Global Maxima Global Maxima
Exceedance Probability
Exceedance Probability
−1 −1
10 10
−2 −2
10 10
p80 − Block Max
−3 −3 p − Block Max
10 10 80
−4 −4
10 10
5 10 15 20 5 10 15 20
OoPBM (MN−m) OoPBM (MN−m)
(a) (b)
Figure 3.6: Short-term distributions for OoPBM for the bin 16 < V < 18 m/s
using 20 simulations for (a) a 20-second block size where independence is verified (see
Table 3.2); and (b) a 5-second size where block maxima are dependent.
Figure 3.6 shows two short-term empirical distributions based on two different
72
block sizes, one of 20-second duration that guarantees independent block maxima and
the other of 5-second duration that almost certainly does not. Short-term distribu-
tions on block maxima for the two block sizes are adjusted based on Eq. 3.6 so as
to be directly represented in terms of exceedance probability in ten minutes. From
Fig. 3.6, it is seen that even with a block length as short as 5 seconds, the tail of
the short-term distribution of OoPBM for the selected 16-18 m/s wind speed bin is
almost identical to the one obtained by using global maxima directly from 20 simu-
lations. This confirms our statement that assumption of independence among block
maxima, even if not justified, has insignificant impact on predicted rare load fractiles.
At the same time, this leads one to the conclusion that empirical distributions based
on the use of larger sized samples of block maxima offer no advantages over the use
of smaller samples of global maxima. It should be pointed out that these findings are
based only on studies with the LE3 loads data sets.
From the preceding discussions, it has been established that adjusting the
block size when extracting extremes data from load time series might sometimes have
little impact on the tails of short-term distributions. Accordingly, now our focus shifts
to whether or not the use of additional simulations, with global maxima extracted
from each simulation, can help to better define distribution tails and, hence, the
73
long-term distribution and consequently the extrapolated rare load. It is of interest
to be able to determine how many global maxima (or simulations) are required to
adequately define rare load fractiles. It is expected that with a larger number of
simulations, additional useful information about rare large loads can be gained. This
should better define short-term load distributions for two reasons. Firstly, additional
simulations will realize more maxima and generally some very large loads that can
help to fill out and better define the tail of the distribution; this will results in
reduced uncertainty in load distributions and, ultimately, in the extrapolated result.
Secondly,with additional simulations, estimation of lower and lower probabilities of
exceedance is possible leading to a reduction in the extent of extrapolation needed
beyond the simulated loads data.
However, the real issue with running simulations is one of practicality. Car-
rying out a large number of simulations can be time-consuming; hence, it would
be beneficial to know what minimum number of simulations is needed to limit the
uncertainty in load predictions to some specified level. A proposal is to enforce a
convergence criterion on the tail fractiles of the empirical short-term distributions,
not on the long-term distributions nor on the extrapolated long-term load itself; for
instance, one could require that the uncertainty in the p-quantile load of the short-
term load distribution must be no larger than some specified amount. This could be
prescribed by requiring a limit on confidence intervals on the p-quantile load. We
propose such a convergence criterion that may be expressed mathematically in the
following manner:
L̂α,p − L̂1−α,p q
< (3.10)
L̂p 100
where the variable, L̂p , in the denominator in Eq. 3.10 represents the empirical es-
timate of the p-quantile load (similar to lp defined in Eq. 3.7 earlier) while the nu-
74
merator represents the (2α − 1)% confidence interval on the p-quantile load. The
right-hand side of the inequality contains the variable, q, which represents the maxi-
mum acceptable percent error permitted on the normalized confidence interval (where
normalization is with respect to the p-quantile load itself). The convergence criterion,
as stated, implies that if the normalized confidence interval based on a certain number
of simulations exceeds the specified maximum acceptable error, additional simulations
will need to be run to reduce this normalized confidence interval. The rationale be-
hind specifying the convergence criterion in this manner is that if the p-quantile is
chosen to be reasonably far in the tail of the short-term distribution, uncertainty in
its estimate can be controlled or limited. Since the tail quite directly influences the
long-term distribution and the extrapolated load, the convergence criterion above,
although it does so only indirectly, aids in the overall purpose of statistical loads
extrapolation. We need to next discuss how the confidence interval in the numerator
of Eq. 3.10 can be estimated from data and how should the maximum percent error,
q, be selected. We need to be cognizant of the excessive level of effort that may be
needed if q is specified very small or if p is specified as a quantile level that is too far
in the tail of the short-term distribution.
Estimates of the confidence interval of the p-quantile load are related to the
number of simulations or data points used in estimating that load. The bootstrap
technique proposed by Efron and Tibshirani [19, 20] makes it possible to estimate
such confidence intervals by a process that involves randomly resampling data, with
replacement, many times. Efron and Tibshirani [19] point out that bootstrapping is a
computer-intensive procedure since the data set may need to be resampled thousands
of times in order to accurately estimate statistics such as confidence intervals on
75
quantities such as quantiles of a distribution. Henderson [32] provides an excellent
review of the bootstrap method and offers many examples of its implementation and
subtleties involved in its use.
Using the bootstrap procedure to form confidence intervals begins with taking
the initial set of data on, say, n global maxima (m1 , m2 , m3 , ..., mn ) and randomly
resampling these data with replacement to form a new set (m1∗ , m∗2 , m∗3 , ..., m∗n ) or a
bootstrap resampling of the same size as the original sample. Note that bootstrap
resamplings will be composed of repeated values from the original sample since, for
each resampling, data are sampled randomly with replacement. The process is re-
peated so as to form a large number, Nb , of bootstrap resamplings. From each of these
sets of n data, individual estimates of the p-quantile can be obtained. From these Nb
estimates, constituting the set, (l1 , l2 , l3 , ..., lNb ), confidence intervals can be found in
the usual manner by ordering the data. These can then be used for the numerator
of Eq. 3.10. The estimate of the p-quantile that is obtained from the original data
represents the denominator of Eq. 3.10.
76
Table 3.4: Estimates of the normalized 90% confidence interval on the 84th percentile
load for OoPBM, FATBM, IPBM and OoPTD based on global maxima from 30
simulations
Wind Speed Normalized 90% Confidence Interval (%)
(m/s) OoPBM FATBM IPBM OoPTD
2<V <4 5.4 8.7 2.2 9.3
4<V <6 6.3 3.2 1.9 6.1
6<V <8 10.8 8.8 3.2 8.8
8 < V < 10 2.4 3.7 4.1 2.5
10 < V < 12 1.8 1.8 1.6 3.8
12 < V < 14 3.2 4.3 2.9 1.4
14 < V < 16 2.6 2.9 2.8 3.2
16 < V < 18 10.5 6.4 1.9 14.3
18 < V < 20 4.6 7.8 6.2 9.9
20 < V < 22 5.2 6.6 4.0 8.9
22 < V < 24 8.5 8.7 6.8 8.7
24 < V < 26 6.1 5.3 5.4 14.0
the benefit of a the larger number is that it will lead to more reliable estimates of
the desired statistic. For the LE3 data set, 5,000 bootstrap resamplings were used to
form confidence interval estimates on load quantiles and these were found to be ade-
quately stable [24]. There is another non-parametric method, based on the binomial
distribution [35], that may be used to estimate confidence intervals. Fogle et al. [24]
showed that the confidence intervals estimates based on the bootstrap and binomial
methods are comparable for the LE3 data set; we will use exclusively the bootstrap
technique in this chapter.
Convergence criteria for four load measures, OoPBM, FATBM, IPBM and
OoPTD were studied using the LE3 loads database. Based on Eq. 3.10, we are in-
terested in computing the percent error in terms of the normalized 90% confidence
77
Table 3.5: Estimates of the normalized 90% confidence intervals on the 84th percentile
global maximum load for four different load types obtained from different numbers
of simulations.
Number OoPBM FATBM IPBM OoPTD
of simulations 16 < V < 18 22 < V < 24 22 < V < 24 16 < V < 18
10 19.0 7.2 2.9 23.6
20 15.0 8.0 5.3 20.7
30 10.5 8.7 6.8 14.3
40 8.9 8.1 6.1 11.6
60 7.2 5.1 3.2 8.4
80 5.6 3.5 2.4 8.9
100 4.2 4.0 2.3 6.3
interval of the 84th percentile ten-minute global maximum for each load type for dif-
ferent numbers of simulations. If the maximum allowable percent error, q, is specified,
the appropriate number of simulations can be run. Table 3.4 shows computed percent
errors when 30 simulations are run for each wind speed bin and for all of the four
loads. The results presented in the table suggest that the convergence criterion is
adequately met if the maximum error permitted is 15% (i.e., if q is equal to 15). For
IPBM and FATBM loads, even a 10% maximum error criterion would be met when
30 simulations are run. For the OoPBM and OoPTD, slowest convergence is seen in
the 16-18 m/s wind speed bin but even there, 30 simulations lead to normalized 90%
confidence intervals on the 84th percentile load that are smaller than 15%.
Short-term load distributions for the four load types are summarized in Fig. 3.7
for the single wind speed bin in each case that showed slowest convergence (or largest
normalized confidence interval) based upon the normalized 90% confidence interval
on the 84th percentile load. The 90% confidence intervals are shown for each empir-
ical distribution at the (1 - 0.84) exceedance probability level. Table 3.5 shows, for
78
0 16 < V < 18 m/s 0 22 < V < 24 m/s
10 10
84th fractile 84th fractile
Exceedance Probability
Exceedance Probability
90% CI 90% CI
−1 −1
10 10
−2 −2
10 10
8 10 12 14 16 18 20 40 60 80 100
OoPBM (MN−m) FATBM (MN−m)
0 22 < V < 24 m/s 0 16 < V < 18 m/s
10 10
84th fractile 84th fractile
Exceedance Probability
Exceedance Probability
90% CI 90% CI
−1 −1
10 10
−2 −2
10 10
6 7 8 9 10 4 6 8 10 12
IPBM (MN−m) OoPTD (m)
Figure 3.7: Short-term distributions for OoPBM, FATBM, IPBM and OoPTD based
on 30 simulations. Bins selected have largest 90% relative confidence bounds on the
84th percentile load as computed based on bootstrap techniques.
controlling wind speed bins, that the normalized confidence interval decreases if more
simulations are run. Therefore, one could choose, if desired in a design situation,
a different level of acceptable error than 15%, and find the appropriate number of
simulations to satisfy that criterion.
79
Table 3.6: Estimates of the normalized 90% confidence intervals on the 84th percentile
global maximum load for four different load types as obtained from block maxima
distributions using different block sizes.
Block OoPBM FATBM IPBM OoPTD
size (sec) 16 < V < 18 16 < V < 18 22 < V < 24 16 < V < 18
5 11.2 8.5 6.8 12.6
10 11.2 8.5 6.8 14.1
15 11.3 8.6 6.8 12.6
20 11.3 8.6 6.8 11.1
30 11.2 8.5 6.8 15.1
60 11.2 8.0 6.8 15.9
block sizes are summarized in Table 3.6; again, we only summarize results for the
single wind speed bin for each load type that resulted in the largest load, as we
did in Table 3.3. It is observed from Table 3.6 that estimates of normalized 90%
confidence intervals on the 84th percentile load are almost the same whether the block
size is selected large, when independence is satisfied, or the block size is small, when
independence is less likely to be satisfied. Moreover, it can be seen that the normalized
confidence intervals for block maxima are almost the same as those for global maxima
(Table 3.4). This observation reaffirms our earlier finding (from Table 3.3) that
extracting more extremes from ten-minute simulations by using block maxima does
not result in significantly new information, compared to simply using global maxima,
at least as far as short-term distributions and extrapolation are concerned.
The preceding discussion suggests that, for the LE3 data set, adequately stable
tails for short-term distributions in all wind speed bins can be obtained if 30 simu-
lations are run. Enforcing the convergence criterion of a maximum percent error of
15% on normalized 90% confidence intervals on the 84th percentile load of short-term
distributions leads us to state that these distributions are reasonably well estimated.
80
3.7 Summary and Conclusions
Using the LE3 data for loads on a model of a 5MW onshore wind turbine, we
have addressed several questions related to the theory and practical implementation
of statistical loads extrapolation. First, we have emphasized the need to understand
which wind speeds tend to cause the largest loads of different types. This information
is useful to have when determining where greater simulation effort is needed.
We have compared the use of global and block maxima, two of several meth-
ods that may be used to extract extremes from a time series, and their differences
in the context of statistical loads extrapolation. Issues related to independence of
block maxima have been addressed using a statistical test of independence. Ignoring
unimportant wind speed bins, it was found that block sizes of around 10-15 seconds
for four loads (OoPBM, FATBM, OoPTD, and IPBM) yielded independent block
maxima when checked at the 1% significant level based on a statistical test. Finally,
though, with the LE3 data, it was demonstrated that there is no advantage gained
from using block maxima over global maxima when short-term loads are estimated.
Moreover, even when block sizes are very small so as to clearly exhibit dependence,
ignoring this again leads to small error in estimation of short-term loads.
81
IPBM) and with 30 simulations, the normalized 90% confidence interval on the 84th
percentile load never exceeded 15%. The rationale behind this objective criterion can
also be used by designers to determine the number of simulations required to obtain
a robust short-term distribution for any chosen limit on the normalized confidence
interval.
82
Chapter 4
4.1 Introduction
The IEC draft design guidelines [43] require statistical extrapolation to predict
long-term loads for offshore wind turbines. Statistical extrapolation, as described
by Eq. 1.1, involves integration of the distribution of turbine loads given specified
environmental states with the likelihood of occurrence of the different environmental
states; the (conditional) load distributions are obtained by means of turbine response
simulations.
Extrapolation methods have been shown to pose several challenges for on-
shore wind turbines [69, 83, 89]; we addressed some issues related to extrapolation
of loads on an onshore wind turbine in Chapter 3. For an offshore wind turbine,
extrapolation methods present several additional challenges [? ]. For one, the off-
shore environment involves, as a minimum, the consideration of waves in addition to
wind; hence, the number of random variables describing the environment increases.
As a result, the domain of integration increases and it can often become impractical
to perform computationally expensive simulations over the entire domain if one uses
the basic extrapolation approach that involves direct integration. It is thus of inter-
est to explore efficient alternative extrapolation techniques for offshore wind turbine
design. A second challenge is that extrapolation of turbine loads needs to recognize
the dependence on two (or more) random processes representing the environment—
83
wind and waves, say—each of which influence turbine loads in distinct ways. Several
studies in recent years have focused on the complexity of these issues in the offshore
environment and have addressed comparisons of alternative methods to extract tur-
bine load extremes [12], possible reduction in simulation effort by careful selection of
critical environmental states [74], the use of the environmental contour method [88],
and the use of a suitable quantile of the wave-related random variable (conditional
on wind speed) in lieu of the full joint wind-wave distribution [98].
84
response statistics for several representative environmental conditions. We then dis-
cuss application of the peak-over-threshold method to derive probability distributions
for turbine loads. We illustrate how long-term loads can be derived using inverse
FORM, first by omitting and then by accounting for turbine load variability (given
environmental state). Comparison of these long-term load predictions based on in-
verse FORM with those obtained by direct integration is discussed. We also discuss
how turbine control actions influence variability in long-term loads and, as a result, a
small number of simulations may not be adequate. Finally, we compare predictions of
rare (long-term) load fractiles based on the POT, the global maxima, and the block
maxima methods.
The primary method for statistical load extrapolation is the direct integration
method (Eq. 1.1), which was discussed in detail in Section 1.3. In this chapter,
we intend to assess the accuracy of more efficient but approximate extrapolation
methods, namely the inverse first-order reliability method and a special case of it, the
environmental contour method. These methods are described below.
In the context of wind turbine loads, suppose our interest is in obtaining es-
timates of the T -year load, lT . For instance, our load of interest might be, say, the
fore-aft tower bending moment at the mudline for an offshore wind turbine. We as-
sume that the uncertainty in the ten-minute global maximum, L, of the selected load
depends on the ten-minute mean wind speed at hub height, V , and the significant
wave height, Hs . From wind turbine response simulations, statistics of L are ana-
85
lyzed and then used to derive the T -year load, lT . Consider a situation where the
joint probability distribution of V , Hs , and L is available in the form of a marginal
distribution for V and conditional distributions for Hs given V , and for L given V
and Hs . Note that if ten-minute time series are simulated to yield statistics on L,
then the T -year load, lT , is expressed in terms of the distribution of L as follows:
If at any point on this sphere, a tangent plane were drawn, the probability of
occurrence of all points on the side of this plane away from the origin is Φ(−β). Since
each point on the sphere is associated with the same reliability index, the desired
load, lT , is also associated with this same reliability. The only difficulty is that the
sphere is in a three-dimensional space of standard normal variables not in the physical
space of the variables, V , Hs , and L. If, however, a mapping from one space to the
other were possible, again, since each point on the sphere is associated with the same
reliability, all the points on the sphere need to be systematically searched until the
largest load is obtained. This largest load is then lT and this procedure is the inverse
first-order reliability method (inverse FORM) [107].
86
The transformation of the random variables involved from the physical X
space (where X = (V, Hs , L)) to the standard normal U space is carried out via the
Rosenblatt transformation [86] such that
where F () denotes the cumulative distribution function in each case. Using this
transformation, one can relate a set of values of the physical random variables, (v, h, l),
to corresponding values of standard normal variables, (u1 , u2 , u3 ), as follows:
¡ ¢ ¡ ¢
u1 = Φ−1 (FV (v)) , u2 = Φ−1 FH|V =v (h) , u3 = Φ−1 FL|V =v,H=h (l) (4.4)
While the above transformation can be done in any order, i.e., u1 , u2 and u3 can be
evaluated independently of each other, the transformation from the standard normal
space back to the physical space, if needed, has to be carried out in the following
order, due to the conditional distributions that are involved:
In other words, v is first computed from the given value of u1 ; this value of v is
then needed to evaluate the inverse conditional distribution of the significant wave
−1
height, FH|V =v , in order to compute h from u2 ; and both these v and h values are
−1
finally required to evaluate the inverse conditional distribution of load, FL|V =v,H=h ,
to compute l from u3 .
The point on the sphere defined by Eq. 4.2 where the load attains its maximum
value is the “design” point, and this load represents the desired T -year return period
load, lT . As an illustration, the sphere in U space and a portion of the corresponding
surface in the physical space are shown in Fig. 4.1; the design point where the load is
87
U2
β
U1
U3
(a) (b)
largest is also identified. Note that this figure is for illustration purposes only, and is
based on arbitrary choices of distribution functions that have no relation to the wind
turbine loads that will be discussed later in this chapter.
88
U2
g(U)=0
gFORM(U)=0
Unsafe region
g(U)<0
Safe region
Design point
g(U)>0
u2
* (u*1,u*2)
E
True limit
state surface
FORM's limit
state surface
u*1 U1
Figure 4.2: True limit state function in the U space, g(U ) = 0, and the correspond-
ing approximate limit state function, gFORM (U ) = 0, according to the first-order
reliability method (FORM). Figure from Saranyasoontorn [87].
finding the design point that will lead to a specified target failure probability, while
FORM can still be used in an iterative way, inverse FORM offers a more efficient
solution. The assumption of a linearized limit state surface is a reason for possible
differences between the long-term load estimated from the inverse FORM and that es-
timated using direct integration (Eq. 1.1); however, this linearization approximation
is not very inaccurate for rare loads associated with very small target probabilities of
exceedance, as has been shown for onshore wind turbine loads [89].
89
U2
U1
(a) (b)
Figure 4.3: Representation of (a) the two-dimensional circle, associated with a re-
liability index, β, in U space, and (b) the associated environmental contour in the
physical space.
This circle, when transformed to the physical space according to Eq. 4.5, is referred to
as the “environmental contour” associated with the target return period, T ; Fig. 4.3
shows an example of such an environmental contour. The design point according to
the EC method is found by searching for the point, on the environmental contour,
where the median load is maximum.
90
The extrapolation methods discussed above require data on load extremes,
which must be obtained from turbine response simulations. Any limitations or ap-
proximations inherent in a simulation model can influence the accuracy of long-term
load predictions. For offshore wind turbines, such model approximation in simula-
tions is introduced in the conventional use of a linear theory to model waves. This
topic will be discussed in detail in Chapter 5. In the present chapter, we focus our
attention on the procedure for extrapolation of offshore wind turbine loads.
A 5MW wind turbine model (Fig. 4.4) developed at NREL [48] closely repre-
senting utility-scale offshore wind turbines being manufactured today is used here. An
onshore version of the same turbine was discussed in Chapter 3 and used in numerical
studies involving statistical load extrapolation. The turbine is a variable-speed, col-
lective pitch-controlled machine with a maximum rotor speed of 12.1 rpm; its rated
wind speed is 11.5 m/s. It is assumed to have a hub height that is 90 meters above
the mean sea level, and a rotor diameter of 126 meters. It is assumed to be sited
91
Figure 4.4: Illustration of the NREL Baseline 5MW offshore wind turbine model
(Picture Credit: Jason Jonkman, NREL).
The linear irregular wave model is the most common model to represent
stochastic ocean waves, and its description can be found in several well-established
92
z
←η(x,t)
MSL z=0 x
← Centerline of
Monopile
z = −h Mudline
Figure 4.5: Schematic illustrating the spatial coordinate, x, the vertical coordinate,
z, and the sea surface elevation, η(x, t), relative to the mean sea level (MSL).
references [3, 11, 90]. We briefly describe this model here. We will limit discussion
to unidirectional (or long-crested) waves. Our interest is in the sea surface elevation,
η(x, t), at a spatial location, x, and at a time instant, t, as illustrated schematically
in Fig. 4.5. Note that the sea surface elevation is measured relative to (and above)
the mean sea level (MSL) where the vertical coordinate, z is zero; also, the centerline
of our turbine monopile will be assumed to be located at x = 0. The sea surface
elevation, according to the linear irregular wave model, is given by the superposition
of linear regular (Airy) waves propagating at different frequencies, as follows.
N
X
η(x, t) = Am cos(ωm t − km x − φm ) (4.7)
m=1
where ωm refers to the frequency of the mth wave component and φm is the associated
random phase assumed uniformly distributed over [0, 2π]. The variable, km , is the
wave number which is related to the frequency, ωm , through the following linear
93
dispersion relation,
2
ωm = gkm tanh(km h) (4.8)
where g is the acceleration due to gravity and h is the water depth. The amplitudes
of the wave components, Am , in Eq. 4.7, are Rayleigh distributed random variables
whose mean square value is given as
where S(ωm ) refers to the one-sided power spectral density function of the sea surface
elevation process. The integer, m, refers to a frequency index that ranges from 1 to
N , the total number of wave components represented in the simulated wave train.
The linear irregular wave model is based on solution of the Laplace equation
in terms of the velocity potential, which is a function of depth, z, as well as of x and
t. The velocity potential, Φ(x, z, t), is given as follows:
N
X cosh(km (h + z))
Φ(x, z, t) = bm sin(ωm t − km x − φm ) (4.10)
m=1
cosh(km h)
The horizontal water particle velocity, u(x, z, t), at depth z (for z ≤ 0), may
be obtained from the spatial derivative of the velocity potential as follows:
∂Φ(x, z, t)
u(x, z, t) = (4.12)
∂x
and the horizontal water particle acceleration, u̇(x, z, t) is obtained by taking the time
derivative of the horizontal water particle velocity.
∂u(x, z, t)
u̇(x, z, t) = (4.13)
∂t
94
Both the horizontal particle velocity and acceleration are required to compute hydro-
dynamic forces on a submerged body, the turbine monopile in this case.
Equation 4.7 is an expression for the simulation of the sea surface elevation,
η(x, t), at any spatial location, x. We are interested in simulating the sea surface
elevation only at the centerline of the monopile, which we arbitrarily choose to locate
at x = 0 (see Fig. 4.5), i.e., we simulate η(0, t); henceforth, we will denote η(0, t)
simply as η(t). Furthermore, for digital simulation of the sea surface elevation process
on a computer in discrete form, we replace continuous time, t, by a discrete time
array, tp = p∆t, where ∆t = T /N , such that p = 1, 2, ..., N . The variable, T , is the
duration of the simulation. We also construct a frequency array, ωm = m∆ω, where
m = 1, 2, ..., N . By virtue of the periodicity of Fourier series (we actually simulate a
duration, T , of a process with period T ) [80], we have
2π
∆ω = (4.14)
T
The sampling theorem [76] requires that the sampling frequency, ωs = 2π/∆T , must
be at least twice the cut-off frequency, ωcut . If ωcut = M ∆ω, then we must satisfy
N ≥ 2M (4.15)
We can now write η(t), according to Eq. 4.7, in terms of discrete time and frequency
as:
N
X
η(tp ) = Am cos(ωm tp − φm ). (4.17)
m=1
95
The summation may further be rewritten as:
" N #
X
η(tp ) = < Am exp(−iφm ) exp (i(m∆ω)(p∆t)) (4.18)
m=1
√
where i = −1 and <[·] denotes the real part of the argument. Recognizing that
∆ω∆t = 2π/N , and defining
we have " N
X µ ¶#
2πm
η(tp ) = < X(ωm ) exp i p (4.20)
m=1
N
Note that the term inside the summation in Eq. 4.20 is the definition of the inverse
discrete Fourier transform of X(ωm ), which can be efficiently computed using inverse
Fast Fourier transform (IFFT) techniques [76], such that
96
where =[·] denotes the imaginary part of the argument, and
Equations 4.22 and 4.23 also imply that the particle velocity and acceleration are
simulated using their respective spectra, Su (ωm ) and Su̇ (ωm ), which are related to
the spectrum, S(ωm ), of the sea surface elevation through their transfer functions
such that
The capability for simulation of linear irregular waves, as outlined above, is imple-
mented in the wind turbine response simulation code, FAST [47], developed at NREL.
We are interested in computing hydrodynamic forces per unit length, f (0, z, t),
at any node located at a depth z, along the centerline of the monopile, which we
arbitrarily chose at x = 0 (Fig. 4.5); for simplicity of notation, we will denote f (0, z, t)
as f . We use Morison’s equation [3] to compute the hydrodynamic forces, f , as follows.
1 πD2
f = fD + fM = CD ρDur |ur | + CM ρu˙r (4.30)
2 4
97
where fD and fM are drag and inertia forces, respectively. The the drag coefficient,
CD , is taken to be 1.0, and the inertia coefficient, CM , is taken to be 2.0 [48]. The
variables, ur and u˙r , are the relative horizontal velocity and acceleration, respectively,
between the water particles and the structure, ρ is the density of water, and D is the
monopile diameter. Note that the water particle velocity and acceleration, computed
using Eqs. 4.22 and 4.23, are valid only up to the mean sea level, and to account for the
effect of the instantaneous water surface, the water particle velocity and acceleration
are modified according to the Wheeler stretching technique [3, 105], which is an
empirical correction.
We are interested in the response of the NREL 5MW offshore wind turbine
while it is in an operating state. Accordingly, we perform response simulations for
mean wind speeds ranging from cut-in to cut-out wind speeds (i.e., 3 m/s to 25 m/s,
here). As a function of the mean wind speed in each simulation, the turbulence in-
tensity is assumed per IEC Class I-B site conditions using the Normal Turbulence
Model [42]. For the hydrodynamic loading on the support structure, irregular lin-
ear waves are simulated using a JONSWAP spectrum [16], employing the procedure
outlined in Section 4.3.2. Hydrodynamic forces are computed using Morison’s equa-
tion; Wheeler stretching [3] is used to account for the influence, above the mean
water level, of the instantaneous wetted surface on kinematics and hydrodynamics.
Stochastic time-domain simulations of the turbine response are performed using the
computer program, FAST. The peak spectral period for the waves is modeled as a
function of significant wave height based on one year’s wave data from the National
Data Buoy Center’s Buoy 44007, where the water depth is 19 meters. We discretize
98
(a) (b)
Figure 4.6: Variation with mean wind speed and significant wave height of the mean of
the maximum values from six simulations of (a) the out-of-plane blade root moment;
and (b) the fore-aft tower base moment.
the domain of the two environmental random variables—the mean wind speed, V ,
and the significant wave height, Hs —using a two-dimensional interval or bin of 2
m/s for mean wind speed and 1 m for significant wave height. We will focus on the
out-of-plane blade moment (OoPBM) at the blade root and the fore-aft tower bend-
ing moment (TBM) at the mudline as the two turbine load variables whose extreme
values are of interest in this study.
99
almost linearly with wave height.
100
Wind Speed (m/s)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(a)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(b)
101
Wind Speed (m/s)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(c)
Figure 4.7: Representative time series of wind speed, sea surface elevation, out-of-
plane bending moment at the blade root, and tower bending moment at mudline for
mean wind speeds of (a) 3.7 m/s, (b) 12.1 m/s and (c) 24.2 m/s. The significant wave
height is fixed at 4.2 m.
102
3
10
V = 24.2 m/s
1
10
0 V = 12.1 m/s
10
−1
10
−2
10
V = 3.7 m/s
−3
10
−4
10
0 0.5 1 1.5 2 2.5
Frequency (Hz)
(a)
4
10
Power Spectral Density ((MN−m)2/Hz)
2 1P
10 Blade 1st Flap.
2P
Blade 1st Edge. nd
Blade 2 Flap.
0
3P
10
V = 24.2 m/s
−2
10
−4
10
0 0.5 1 1.5 2 2.5
Frequency (Hz)
(b)
103
4
10
Tower 1st
0
10
−2
10 V = 12.1 m/s
V = 3.7 m/s
−4
10
0 0.5 1 1.5 2 2.5
Frequency (Hz)
(c)
Figure 4.8: Variation of power spectral density with mean wind speed for (a) wind
speed, (b) out-of-plane bending moment at the blade root, and (c) tower bending
moment at mudline. The significant wave height is fixed at 4.2 m.
104
Wind Speed (m/s)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(a)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(b)
105
Wind Speed (m/s)
20
0
0 50 100 150 200
Wave Elevation (m)
5
0
−5
0 50 100 150 200
Out−of−plane Blade Root Momemt (MN−m)
15
5
−5
0 50 100 150 200
Tower Base Momemt (MN−m)
150
50
0
0 50 100 150 200
time (sec)
(c)
Figure 4.9: Representative time series of wind speed, sea surface elevation, out-of-
plane bending moment at the blade root, and tower bending moment at mudline for
significant wave heights of (a) 0.5 m, (b) 4.2 m and (c) 9.4 m. The mean wind speed
is fixed at 12.1 m/s.
106
4
10
Hs = 9.4 m
st
Tower 1
−2
10
−4
10
0 0.5 1 1.5 2 2.5
Frequency (Hz)
Figure 4.10: Variation of power spectral density of fore-aft tower bending moment
with significant wave height. The mean wind speed is fixed at 12.1 m/s.
The effect of waves is studied by comparing the turbine response for signif-
icant wave heights of 0.5 m, 4.2 m and 9.4 m, while the mean wind speed is held
constant at 12.1 m/s. Figure 4.9 clearly shows larger peaks in the TBM time series
with increasing wave height. Blade loads are seen to be insensitive to wave height
variation. Accordingly, power spectra for tower loads alone are presented in Fig. 4.10
where no significant influence of change in wave height is noted, except at frequencies
below 0.2 Hz where wave energy is dominant, with sea surface elevation peak spectral
frequencies being 0.14, 0.10 and 0.08 Hz for significant wave heights of 0.5, 4.2 and
9.4 m, respectively.
107
Table 4.1: Ten-minute statistics of out-of-plane bending
moment at the blade root for different wind speed and
wave height bins.
Out-of-plane moment at blade root
V Hs
Max Mean SD Skew. Kurt. PF
(m/s) (m)
(all in MN-m) -
12.0 0.5 12.5 8.2 1.5 -0.05 2.64 2.85
12.0 4.5 12.7 8.3 1.6 -0.25 2.74 2.80
12.0 9.5 12.2 8.2 1.5 -0.14 2.61 2.60
4.0 4.5 4.5 2.3 0.6 0.49 2.83 3.65
12.0 4.5 12.7 8.3 1.6 -0.25 2.74 2.80
24.0 4.5 9.3 3.0 2.0 -0.14 2.89 3.07
Note: V : Mean wind speed, Hs : Significant wave height,
Max: Ten-minute maximum, SD: Standard deviation,
Skew.: Skewness; Kurt.: Kurtosis, PF: Peak factor =
(Max - Mean)/SD.
108
Tables 4.1 and 4.2 summarize statistics of the blade and tower loads, respec-
tively, obtained from six simulations each for the wind and wave conditions discussed
in Figs. 4.7-4.10. OoPBM statistics, as was discussed before, are insensitive to wave
height. Also, the mean and maximum OoPBM are systematically higher around
the rated wind speed compared to other wind speeds, e.g., the maximum OoPBM
load near the rated wind speed is about 35% larger than that near the cut-out wind
speed. The standard deviation of the OoPBM near the rated wind speed is, however,
about 20% smaller than that near the cut-out wind speed. Tables 4.1 and 4.2 also
present values of skewness, kurtosis and peak factors, which give insights into the
non-Gaussian character of the load processes. Assuming that a process is Gaussian,
we can estimate a theoretical peak factor from the number of local maxima in a
realization—say, for example, a ten-minute simulation—of the random process [62].
For the OoPBM process, the value of the Gaussian peak factor for all three wind
speeds and wave heights is about 3.3, which is larger than the computed peak factors
(shown in Tables 4.1 and 4.2) when the skewness is negative, and smaller when the
skewness is positive (only for the case where V = 4 m/s and Hs = 4.5 m). Note
that the kurtosis values for the OoPBM process are also different from those for a
Gaussian process, whose kurtosis value is 3.0. The implication of the non-Gaussian
character of turbine load processes is that the extremes for rare (small) probability
levels might be very different from those predicted by a Gaussian process with the
same mean and variance.
The TBM process statistics are also interesting—when studied with variation
in wave height, it is seen that mean levels change only very slightly. This is because
the sea surface elevation process has zero mean, and changes in significant wave height
do not affect the mean response, especially since no current is assumed present. On
109
the other hand, as the significant wave height is increased from 4.5 m to 9.5 m, the
standard deviation and maximum value of the TBM increase by about 26% and 14%,
respectively. This is because significant wave height, which is a measure of the energy
(variance) of the wave process, directly affects the energy of the turbine response
process, and thus the maximum and variance of the TBM process increase with wave
height. As was the case for the blade loads, peak factors for the TBM process are
different from those for a Gaussian process, and so are the skewness and kurtosis
values. Again, these reflect the non-Gaussian nature of the tower load process.
From the preceding results, it can be concluded that: (1) blade and tower
loads are largest around the rated wind speed, but peak factors are lowest there; (2)
blade loads are independent of wave height; and (3) maximum tower loads increase
systematically with wave height. Hence, turbine long-term extreme loads are expected
to be governed either by mean wind speeds near rated (as the mean and maximum
turbine response are higher there) or by higher-than-rated wind speeds where larger
variability in loads and associated large peak factors could lead to large extreme
values. Also, tower long-term extreme loads are likely to result from larger wave
heights.
The short-term distribution of turbine loads, FL|X (l), which enables predic-
tion of long-term loads according to the direct integration method (Eq. 1.1) or in-
verse FORM (Section 4.2.1), requires data on load extremes. We use the peak-over-
threshold (POT) method here, which is discussed in the IEC guidelines [43], to extract
load extremes from time series data. While the POT method, compared to the use of
global maxima, can provide a larger number of load extremes from a given number of
110
simulations, the load fractile for POT data, equivalent to that for a specified fractile
for the global maxima, is also transformed farther out into the tail of the POT dis-
tribution. Therefore, a relevant question to ask is whether the additional extremes
provided by the POT method are actually advantageous for extrapolation or not;
this question, which also arises with the block maxima method as was discussed in
Chapter 3, is our focus here. In addition, we also discuss how the POT method can
be used with inverse FORM to derive long-term loads for design.
£ ¤n
FL|X =x (l) = FLP OT |X =x (l) (4.31)
where n is the expected number of peaks above the chosen threshold in ten minutes.
Note that this relation is similar to that used for block maxima (Eq. 3.6). Equa-
tion 4.31 is based on the assumption that the peaks above the chosen threshold (for a
given X) are independent. If the load associated with a non-exceedance probability
in ten minutes equal to p is of interest, that load, lp , based on the POT distribution,
is associated with a non-exceedance probability, p1/n , and may be estimated as:
£ ¤
lp = FLP OT |X p1/n (4.32)
111
Note that as the selected threshold level is increased, the expected number of
peaks that are retained, in a given duration (which is ten minute here), decreases, and
when this number is unity, the POT method becomes more and more like the global
maximum method since then, on average, one peak is extracted from each simulation.
For typical threshold levels, the expected number of retained peaks is significantly
larger than unity and, as a result, p1/n can approach values that are close to unity. As
an example, if the expected number of peaks above a chosen threshold is 80, then the
non-exceedance fractile level for POT data corresponding to the ten-minute median
extreme load is 0.51/80 = 0.99137. If loads corresponding to this probability level are
to be established non-parametrically from simulations, at least 1/(1 − 0.99137) or 116
peak values above the chosen threshold must be available for the wind speed and wave
height values represented by X. For tight confidence intervals on estimates of such
rare load levels, the number of peaks (and, thus, the number of simulations) might
even need to be a few orders of magnitude higher. Note that extrapolation may often
be required then for two reasons: (1) to estimate rarer fractiles of the ten-minute
extreme load instead of the median, as the minimum number of required data may
far exceed the amount of POT data that are available from limited simulations; and
(2) to have tight confidence intervals on any predicted POT load fractiles.
112
loads change as the number of simulations is increased, and whether six simulations
are adequate or not, for the offshore wind turbine studied here.
We first estimate long-term loads using the 2-D environmental contour method.
Then, we compare these 2-D loads with those obtained using a full 3-D inverse reli-
ability approach. We start by using six ten-minute turbine response simulations for
each environmental state to establish turbine load statistics, and subsequently inves-
tigate the effect of varying the number of simulations on long-term load predictions.
All the long-term loads discussed hereinafter correspond to a return period of 20 years
which, according to the design standards [42, 43], is the target service life for which
wind turbines are typically designed.
113
I-B wind regime (for which our turbine model is being considered), we assume that
the ten-minute mean wind speed, V , at hub height has an average value of 10 m/s and
that it can be described by a Rayleigh distribution. We truncate this distribution
below 4 m/s and above 24 m/s, since we are interested only in studying turbine
loads during operation. The significant wave height, Hs , conditional on the mean
wind speed, is assumed to be represented by a two-parameter Weibull distribution.
The expected value of Hs given V is based upon the JONSWAP correlation between
wind and waves [11], while a coefficient of variation for Hs given V is assumed to be
constant at 0.2.
9
Significant Wave Height (m)
8
7
6
5
4
3
2
1
0
5 10 15 20 25
Mean Wind Speed (m/s)
114
Table 4.3: Comparison of 20-year blade and tower loads estimated by differ-
ent methods, when load extremes data are based on the peak-over-threshold
approach.
20-year design for OoPBM 20-year design for TBM
Method V Hs OoPBM V Hs TBM
(m/s) (m) (MN-m) (m/s) (m) (MN-m)
EC 12.0 6.2 12.8 12.0 6.2 105.2
EC with corrections 12.0 6.2 13.2 12.0 6.2 107.9
3D Inverse FORM 14.0 5.5 13.6 16.0 5.5 119.9
Note: EC: Environmental Contour, FORM: First-order Reliability Method,
OoPBM: Out-of-plane bending moment at the blade root, TBM: Fore-aft tower
bending moment at the mudline.
mean of the significant wave height conditional on the mean wind speed, larger wave
heights are associated with larger wind speeds in the upper right-hand corner of the
environmental contour in Fig. 4.11. To estimate the long-term load of interest, we
seek the maximum value of the median ten-minute extreme turbine load given X on
this environmental contour. This median load is mapped to a corresponding fractile
based on POT data by setting p to be 0.5 in Eq. 4.32. Estimated 20-year loads are
presented in Table 4.3. The 20-year OoPBM load is 12.8 MN-m which is associated
with a mean wind speed of 12 m/s and a significant wave height of 6.2 m. The 20-
year TBM load is 105.2 MN-m and it also results from the same wind speed and wave
height. That the “design” wind speed is close to the rated wind speed is expected
as median extreme turbine loads are largest there, as was discussed earlier (Fig. 4.6).
The design wave height of 6.2 meters is the larger of the two possible wave heights
on the 20-year environmental contour that accompany a mean wind speed of 12 m/s.
115
Table 4.4: Required fractiles for the design environmental states with the 2-D envi-
ronmental contour method.
Average Required Largest
Load number of fractile empirical
peaks, n 0.51/n fractile
OoPBM 87.2 0.9921 0.9981
TBM 80.2 0.9914 0.9970
tion, given the number of peaks above the threshold retained from the six simulations.
Table 4.4 shows that for both blade and tower loads, the required high non-exceedance
POT fractiles are smaller than the highest available empirical non-exceedance frac-
tile, 1 − 1/(6n + 1). This suggests that extrapolation is not necessary to arrive at the
turbine long-term loads with the EC method; nevertheless, the method has accuracy
limitations both because it does not employ the full distribution of turbine loads con-
ditional on X and because even a non-extrapolated fractile is subject to statistical
uncertainty due to limited data. A separate source of inaccuracy, not related to the
quantity of data used, is due to the assumed linear failure surface in any first-order
reliability method procedure such as the EC method. To assess this latter inaccuracy,
we estimate long-term loads using direct integration. We model the conditional load
distribution in Eq. 1.1 as a step function that attains a unit value at the median
load, as discussed in Section 4.2.2. To yield the desired target probability of fail-
ure, the long-term loads are found to be 12.7 MN-m and 110.2 MN-m, respectively,
for OoPBM and TBM, which are very close to those obtained by the environmental
contour method. Hence, we conclude that the EC method is not grossly inaccurate
relative to an exact integration approach that works with the same data. Still, there
are other reasons why the EC-based long-term loads may not be correct; these reasons
have to do with incomplete modeling of the conditional distribution of turbine loads
116
given wind speed and wave height. This is addressed next.
While the full distribution of loads (given environmental conditions) can be em-
ployed in a 3-D inverse FORM approach, this requires far more computational effort.
An alternative strategy is to apply a correction to the 2-D EC long-term loads [106],
as has been demonstrated for an onshore wind turbine [89]. In this approach, the
neglected response variability in the EC method is attributed to (1) background vari-
ability in the median extreme response, L̂, which accounts for the variability in the
median response with changing environmental states (corresponding to different re-
turn periods); and (2) response variability which arises due to different stochastic
components for a given environmental condition (this variability is modeled by quan-
tifying different load fractiles at the same environmental conditions). A corrected
short-term extreme response is expressed as L = L̂², where ² is a unit-median ran-
dom variable and it represents the variability in the extreme response at a given set
of environmental conditions. Standard deviations of the natural logarithms of these
two random variables, L̂ and ², are given by
³ ´
ln L̂T1 /L̂T2 ln (²p1 /²p2 )
σln L̂ = ; σln ² = −1 . (4.33)
βT1 − βT2 Φ (p2 ) − Φ−1 (p1 )
where T1 is the target return period (20 years); T2 is a slightly smaller return period
(taken as 16 years here); and βT1 and βT2 are associated reliability indices. Also, p1
is the median fractile level and p2 is a higher fractile (taken as the 0.83 fractile here)
fractile level, while ²p1 and ²p2 are associated load values from the simulations at the
selected environmental conditions associated with the T -year long-term load. The
117
correction factor, R, is expressed as:
LT
R= = exp [(σln L − σln L̂ ) βT ] (4.34)
L̂T
2 2 2
where σln L = σln L̂ + σln ² .
After applying this correction, long-term loads are estimated to be 13.2 MN-
m and 107.9 MN-m for OoPBM and TBM, respectively (see Table 4.3). Response
variability is largely responsible for the change in long-term loads here. However,
these corrected blade and tower long-term loads are only about 3% larger than the
corresponding 2-D EC values.
If, instead of seeking only the median extreme load given X, the full prob-
ability distribution of the turbine extreme load, L, is established by simulations, a
search is needed for the maximum value of a different fractile, p3 , on load extremes
consistent with each environmental state (V, Hs ) and with the specified target return
period, T , or target probability, PT (and, hence, the associated reliability index, β,
where PT = Φ(−β)). The desired fractile level is p3 = Φ(u3 ), and by rewriting Eq. 4.2
in terms of u3 , and expressing u1 and u2 in terms of the physical variables using the
Rosenblatt transformation (Eq. 4.4), we have the following expression for p3 :
µq ¶
2 £ ¡ ¢¤2
p3 = Φ β 2 − [Φ−1 (FV (v))] − Φ−1 FH|V (h) (4.35)
where Φ() and Φ−1 () refer to the standard normal cumulative and inverse cumulative
distribution functions, respectively, and F () denotes the respective cumulative dis-
tribution function of environmental random variables. For the POT data, the load
fractiles are estimated according to Eq. 4.32.
118
With this 3-D inverse FORM approach, the 20-year long-term loads (here ob-
tained by searching only on gridded (V, Hs ) pairs where simulations were carried out)
are 13.6 MN-m and 119.9 MN-m for the blade bending moment and tower bending
moment, respectively (see Table 4.3). These same loads are obtained using the direct
integration method, which establishes the accuracy of the 3-D inverse FORM results.
These 3-D (inverse FORM) long-term loads are roughly 6% and 14% larger, respec-
tively, for the blade and tower bending moments than those obtained with the 2-D
EC method. Interestingly, the design wind speed with the 3-D formulation, which is
14 m/s for the blade load and 16 m/s for the tower load, is no longer near the rated
wind speed, as was the case with the 2-D formulation. This suggests that accounting
for the full conditional load variability (as a function of wind speed and wave height)
is important.
Note that for a pitch-controlled turbine, the type of turbine considered here,
the rated wind speed is sometimes assumed to be the design wind speed that controls
long-term loads. Furthermore, in order to reduce the simulation effort, a limited
number of wind speeds (as few as three), close to rated, may sometimes be selected
for simulations, as suggested in Annex G of the draft IEC guidelines for offshore wind
turbines [43]. If we employed such an approach, with a small wind speed bin size of 1
m/s, we would not include the design wind speeds of 16 m/s in our simulations that
were found here to control the 20-year load.
We should note that the long-term loads based on the 2-D EC method and on
the 3-D inverse FORM procedure were calculated using simulations for a discrete set of
gridded values of V and Hs . In subsequent discussions, we focus on the environmental
state (i.e., the V and Hs values) at the 3-D design point in Table 4.3. For the blade
and tower loads, these environmental states corresponding, respectively, to V = 14
119
Table 4.5: Required fractiles for design environmental states
with the 3-D inverse FORM approach.
Required Average Required Largest
Load load fractile, number of fractile for empirical
1/n
p3 peaks, n POT, p3 fractile
OoPBM 0.99998997 66.2 0.99999985 0.9975
TBM 0.99999613 74.2 0.99999995 0.9978
Note: OoPBM: Out-of-plane bending moment at the blade
root, TBM: Fore-aft tower bending moment at the mudline.
m/s, Hs = 5.5 m and V = 16 m/s, Hs = 5.5 m in the 3-D approach, are studied in
greater detail in the following.
In order to assess the accuracy of estimated p3 -fractile loads for the design
environmental state in the 3-D inverse FORM approach, it is useful to first determine
if the required fractile needed to be extrapolated from the available POT data. Ta-
ble 4.5 shows that for both loads, the required POT fractiles are significantly higher
than the largest available empirical fractile. In our non-parametric model for the
load distributions, we assume saturation of the tail and somewhat simplistically es-
timate the required fractile as the largest observed value in the simulations. The
load distributions (Fig. 4.12), however, show that the assumption of saturation of
tails (by saturation, we mean that the load level for smaller exceedance probability
levels would be the same as the largest observed load value) cannot be justified. As
a result, 3-D inverse FORM long-term loads based on this non-parametric approach
here would clearly be low and unconservative. We will seek a more accurate estimate
of these long-term loads by establishing short-term load distributions, at the design
environmental states, with more stable tails; this will in turn require that we run
more simulations.
120
0 0
10 10
Exceedance Probability
Exceedance Probability
−2 Related to Median −2 Related to Median
10 10
in 10−minutes
in 10−minutes
−4 −4
10 10
−6 −6
10 10
Related to p3 fractile
Related to p3 fractile
−8 −8
10 10
0 5 10 15 20 0 50 100 150 200
Out−of−plane Blade Root Moment (MN−m) Tower Base Moment (MN−m)
(a) (b)
Figure 4.12: Load distributions of the POT data for governing environmental states
based on 6 simulations for (a) a mean wind speed of 14 m/s and significant wave
height of 5.5 m for OoPBM; and (b) a mean wind speed of 16 m/s and a significant
wave height of 5.5 m for TBM.
An interesting aspect that may be seen from Fig. 4.12 is that the maximum
observed load value, which has a strong bearing on the extrapolated load, is signifi-
cantly larger than the other load values in the tail of the distribution. To investigate
what conditions bring about this large load, we study in Fig. 4.13 relevant time histo-
ries of wind speed, OoPBM, and TBM time series along with that of the blade pitch
angle for a single ten-minute simulation (out of six for the same V − Hs pair) that
included this large load. The maximum loads are seen to occur when the blade pitch
angle suddenly reduced to zero at time instants of roughly 20 sec, 100 sec, and 175
sec. This is due to the control system for this pitch-controlled turbine, which is such
that the blades start to pitch when the instantaneous wind speed exceeds the rated
wind speed of 11.5 m/s. At instants when the wind speed falls below rated, the pitch
angle reduces to zero, and if the wind speed increases before the blade can pitch back,
large loads result.
121
V13H5N9
Wind Speed (m/s)
16
12
8
0 50 100 150 200
Blade Pitch (degrees)
10
5
0
0 50 100 150 200
Out−of−plane Blade Root Moment (MN−m)
14
10
6
2
0 50 100 150 200
Tower Base Moment (MN−m)
110
80
50
20
0 50 100 150 200
time (sec)
Figure 4.13: Time series of wind speed, blade pitch, OoPBM, and TBM for a mean
wind speed of 14 m/s and a significant wave heights of 5.5 m.
Since these large loads due to control actions are observed only in one out of six
simulations, the distribution tails may only saturate and have better definition than in
Fig. 4.12 if more such large load values result upon performing additional simulations.
We therefore carry out additional simulations for the governing environmental states
and find that at least 60 and 150 simulations, respectively, are needed for the blade
and tower loads. The corresponding distributions, shown in Fig. 4.14, also illustrate
how the distribution tails fill in and hence become more reliable. Clearly, due to blade
pitch-control actions, performing only six simulations per environmental state may
be inadequate to obtain reliable distributions by means of parametric model fits to
the data; this is why non-parametric fractiles were employed with the 2-D and 3-D
122
First 6 Simulations First 6 Simulations
Next 54 Simulations Next 144 Simulations
0 0
10 10
Exceedance Probability
Exceedance Probability
−2 Related to Median −2 Related to Median
10 10
in 10−minutes
in 10−minutes
−4 −4
10 10
−6 −6
10 10
Related to p3 fractile
Related to p3 fractile
−8 −8
10 10
0 5 10 15 20 0 50 100 150 200
Out−of−plane Blade Root Moment (MN−m) Tower Base Moment (MN−m)
(a) (b)
Figure 4.14: Load distributions of the POT data for governing environmental states
based on (a) 60 simulations of OoPBM for a mean wind speed of 14 m/s and a
significant wave height of 5.5 m; and (b) 150 simulations of TBM for a mean wind
speed of 16 m/s and a significant wave height of 5.5 m.
With the more reliable POT load distribution tails now possible as a result of
the larger number of simulations, we attempt fits with parametric models. With a
two-parameter Weibull distribution on the tails and a least squares basis for fitting,
Figs. 4.15(a) and 4.15(b) show that for the required fractiles of Table 4.5, loads of 15.3
MN-m and 147.1 MN-m result for the blade and tower loads, respectively. These loads
are about 13% and 23% larger for blade and tower bending moments, respectively,
than those from the non-parametric approach and based on only six simulations. This
is expected due to the unconservatively assumed saturation of the non-parametric
distribution tails before. As seen, a large number of simulations are required to yield
reliable distribution tails and more accurate long-term loads. Note that the accuracy
of the turbine response simulations and, hence, that of the load distributions and
predicted long-term loads also depends on the model uncertainties associated with
123
1e−008 1e−008
Related to p3 fractile Related to p3 fractile
1e−005 1e−005
0.001 0.001
0.01 0.01
0.1 0.1
0.3 0.3
0.5 0.5
0.7 0.7
Unused Data Unused Data
Used Data Used Data
Weibull Fit Weibull Fit
0.9 0.9
8 12 16 20 50 100 150 200
Out−of Plane Blade Root Moment "M" (MN−m) Tower Base Moment "M" (MN−m)
(a) (b)
Figure 4.15: Two-parameter Weibull distribution fits to POT data for (a) OoPBM
(for V = 14 m/s, Hs = 5.5 m), and based on 60 simulations (b) TBM (for V = 16
m/s, Hs = 5.5 m) and based on 150 simulations.
the simulation code (FAST). We, however, do not focus on such uncertainties of the
simulation software nor on the turbine model used; as such, this is a limitation of
this study. In Chapter 5, we will address inaccuracies that result due to the use of
simplistic (linear) wave models that may not be appropriate at shallow water depths.
124
1e−008 1e−008
1e−005 1e−005 p3 fractile
p3 fractile
0.001 0.001
0.01 0.01
0.1 0.1
0.3 0.3
0.5 0.5
0.7 0.7
Unused Data Unused Data
Used Data Used Data
Weibull Fit Weibull Fit
0.9 0.9
8 12 16 20 50 100 150 200
Out−of Plane Blade Root Moment "M" (MN−m) Tower Base Moment "M" (MN−m)
(a) (b)
Figure 4.16: Two-parameter Weibull distribution fits to global maxima data for (a)
OoPBM (for V = 14 m/s, Hs = 5.5 m), and based on 60 simulations (b) TBM (for
V = 16 m/s, Hs = 5.5 m) and based on 150 simulations.
We extract global maxima, the maximum values in each ten-minute turbine re-
sponse time series for the design environmental states, and fit two-parameter Weibull
distributions to the tails of these global maxima data. From these distributions,
we then estimate the load fractiles required with the 3-D inverse FORM approach.
Figures 4.16(a) and 4.16(b) show these Weibull fits for the blade and tower loads,
respectively. Long-term (20-year) load predictions obtained using the global maxima
method are 14.5 MN-m and 136.6 MN-m for the blade and tower bending moments,
respectively. These are only about 5% and 7% smaller for the blade and tower loads,
respectively, than those obtained with the POT method with parametric distribution
fits (Fig. 4.15). The slightly larger differences in 20-year load predictions for the
tower bending moments are likely due to relatively poor distribution fits with both
methods.
An important issue when using the POT method is related to the choice of
threshold. As the threshold level is increased, the number of peaks decreases and they
125
Table 4.6: Effect of threshold level on the estimate of load fractiles for
the design environmental state for OoPBM (V = 14 m/s, Hs = 5.5 m,
where p3 = 0.99998997) and TBM (V = 16 m/s, Hs = 5.5 m, where p3 =
0.99999613). Threshold = Mean + Nσ ×(Standard deviation).
OoPBM TBM
Threshold
Exceedance Load Exceedance Load
level, Nσ n probability fractile n probability fractile
1/n 1/n
1 − p3 (MN-m) 1 − p3 (MN-m)
1.4 66.8 1.50×10−7 15.3 85.6 4.52×10−8 147.1
1.7 44.7 2.24×10−7 14.9 56.4 6.86×10−8 148.1
2.0 28.9 3.47×10−7 14.8 34.9 1.11×10−7 145.9
2.3 17.1 5.86×10−7 14.7 20.2 1.91×10−7 143.0
2.7 8.1 1.23×10−6 14.6 9.7 3.97×10−7 139.7
3.0 4.0 2.51×10−6 14.6 4.9 7.97×10−7 142.4
MAX 1.0 1.00×10−5 14.5 1.0 3.87×10−6 136.6
Note: n: Average number of peaks in ten minutes, OoPBM: Out-of-plane
bending moment at the blade root, TBM: Fore-aft tower bending moment
at the mudline.
become more independent1 . For an appropriately high threshold, the POT method
may result in the same number of load extremes, on average, as the global maxima
method. We now estimate required fractiles for the governing environmental state
with the 3-D inverse FORM approach using different thresholds and two-parameter
Weibull parametric fits to the POT distribution tails. Table 4.6 shows computed
fractiles for blade and tower loads. For the blade loads, it can be seen that the
variation in load fractiles with different thresholds is not significant. The reason is
that very good parametric fits, such as those shown in Fig. 4.15(a), are obtained
for all threshold levels. For the tower loads, the required load fractiles show slightly
greater variation with different threshold choices which may be partly due to less
1
Note that Blum’s test for independence (Section 3.5.1), which was used for block maxima, cannot
be used for POT as extremes with the POT method tend to occur in clusters.
126
Table 4.7: Average plus one standard deviation values of B from block
maxima with different block sizes based on 60 simulations for the out-
of-plane bending moment at the blade root (OoPBM) and based on 150
simulations for the fore-aft tower bending moment at mudline (TBM).
B values
Block size OoPBM TBM
(sec) (V = 14 m/s, Hs = 5.5 m) (V = 16 m/s, Hs = 5.5 m)
5 28.39 6.68
10 12.84 5.77
15 8.34 3.89
20 5.62 3.67
25 4.17 3.42
30 3.96 3.18
60 2.54 2.82
Note: The critical value, Bcr , at the 1% significance level is 4.23.
evident stable trends in the distribution tails for these loads (as seen, for example, in
Fig. 4.15(b)). Based on these observations, we conclude that the agreement between
predicted long-term loads using the POT and global maxima methods is generally
good and is independent of threshold choice as long as the distribution tails are well-
defined enough to allow a good parametric fit.
Block maxima, as discussed in Section 3.5, are alternative load extreme statis-
tics that may be extracted from a time series. This is done by partitioning the time
series into non-overlapping blocks of constant duration and from each of these blocks,
the largest value is extracted as the extreme. The block size directly affects whether
the extracted extremes are independent or not, and we presented an objective crite-
rion, based on the B statistics (Eq. 3.9) according to Blum’s test, to establish the
independence of extremes and for selection of an optimal block size. To obtain a sta-
127
ble estimate, the average of B values over multiple simulations (or the average plus
one standard deviation value, which is a stricter criterion) was used. The average plus
one standard deviation values of the B statistic for the design environmental states
for blade and tower loads are shown in Table 4.7. It is found that at the 1% signifi-
cance level, when the critical B value is 4.23, block sizes of 25 seconds and 15 seconds
would satisfy the independence criterion for blade and tower loads, respectively. As
a conservative measure we choose the same block size of 30 seconds, which results in
20 extremes per ten-minute time series, for both the blade and tower loads. To esti-
mate the 20-year load according to the inverse FORM approach, we fit two-parameter
Weibull distributions to the tails of the block maxima data, and compute the required
fractile. Long-term loads obtained using the block maxima data are 15.1 MN-m and
143.2 MN-m for the blade and tower loads, respectively. These predictions are very
close to those obtained using the POT and global maxima approaches (Table 4.6).
This observation is in line with our earlier conclusion in Chapter 3, where we showed
that rare fractiles predicted by the block maxima method, irrespective of the block
size, are very close to those predicted by the global maxima method. Therefore, it
can be concluded that long-term loads predicted using the three extreme models—
POT, global maxima and block maxima—are generally similar if a sufficient number
of simulations are available to obtain reliable distribution tails.
128
Table 4.8: Estimates of normalized 90% confidence interval on the 84th per-
centile load for design environmental states for the out-of-plane bending mo-
ment at the blade root (OoPBM) and the fore-aft tower bending moment at
the mudline (TBM).
OoPBM (V = 14 m/s, Hs = 5.5 m) TBM (V = 16 m/s, Hs = 5.5 m)
Normalized 90% CI Normalized 90% CI
Number of (in %) Number of (in %)
simulations Global Block POT simulations Global Block POT
6 8.0 7.9 7.7 6 18.8 21.0 21.4
12 9.4 9.2 8.6 20 21.6 21.8 21.7
20 5.4 4.0 3.7 50 6.8 5.7 5.7
30 3.2 3.0 3.0 60 5.4 5.1 5.1
40 3.1 3.2 3.0 80 4.7 3.9 3.9
50 2.9 2.9 2.9 100 4.5 3.9 3.9
60 3.0 2.4 2.4 120 3.9 3.9 3.8
Note: CI: Confidence interval. Normalized CI = (L95,84 − L5,84 )/L84 , where
L84 is the 84th percentile; L5,84 and L95,84 are the 5th and the 95th percentiles
on estimates of the 84th percentile obtained from 1,000 bootstrap resamplings.
confidence interval (see Eq. 3.10), defined as the 90% confidence interval divided by
the point estimate of the 84th fractile, is below an acceptable limit, we claim that the
number of simulations are enough to obtain a reliable distribution. Such normalized
90% confidence intervals for the governing environmental states for blade and tower
loads, obtained from 1,000 bootstrap resamplings, are summarized in Table 4.8. It
is observed that for a given number of simulations, the normalized confidence inter-
vals are very close for the global maxima, block maxima, and POT methods; this is
because the fractile levels for the global maxima and the other methods are related
(Eqs. 3.7 and 4.32). The confidence interval narrows as the number of simulations
is increased. Moreover, it takes a larger number of simulations for the tower load,
compared to that for the blade load, to obtain any specified small confidence interval,
a fact that is consistent with our earlier discussion (see Fig. 4.14) where we showed
129
that to obtain stable tails, more simulations are required for the tower loads. If we
choose a normalized confidence interval of 5% to be acceptable, then at least 30 and
80 simulations are required for the blade and the tower loads, respectively. We chose
to use an even larger number of simulations (60 and 150 for blade and tower loads,
respectively) to obtain better defined distribution tails and more robust predictions
of turbine long-term loads.
Our objective in this study was to derive long-term loads for a utility-scale
5MW offshore wind turbine sited in 20 meters of water. The focus was on the out-of-
plane blade bending moment at a blade root and the fore-aft tower bending moment
at the mudline. Load extremes data needed to establish short-term load distributions
were extracted from time series of turbine response simulations using the peak-over-
threshold, block maxima and global methods. Long-term loads were estimated using
2-D and 3-D inverse first-order reliability method approaches (the former is also re-
ferred to as the environmental contour or EC method) and were compared with direct
integration. The following are general conclusions for the offshore wind turbine stud-
ied:
130
than those based on higher-than-median fractiles (3-D inverse FORM). The 3-D
inverse FORM is found to be as accurate as the direct integration method and
is significantly more efficient.
• Long-term load predictions based on the global maxima, block maxima and
POT methods were close as long as distribution tails were reliable and well
defined.
• Importantly, the design wind speed governing turbine long-term loads is not the
rated wind speed (as is often the case) but is above rated when load variability
is considered.
• A chief source of load variability results from blade pitch control actions that
result in large loads such that the tails of the short-term load distribution are
not reliable unless a large number of simulations is performed.
These conclusions, while particular to the turbine model studied, are useful
to consider for any simulation-based exercise that seeks to predict turbine long-term
loads. This study also suggests that the effect of control actions on extreme loads
needs careful study; in particular it is of interest to investigate methods to account
for variability due to control actions since they can alter the tails of load distributions
in different ways than for loads that result from uncontrolled turbine states.
131
wind turbine sites. In the next chapter, we investigate the influence of the use of an
alternative nonlinear wave model in establishing turbine long-term loads.
132
Chapter 5
5.1 Introduction
While addressing different load cases according to the IEC guidelines for off-
shore wind turbines [43], designers are required to estimate long-term extreme and
fatigue loads; this is usually done by carrying out time-domain stochastic turbine
response simulations. The simulation of offshore wind turbine response involves, as
discussed in Chapter 4, simulations of the stochastic inflow wind field on the rotor
plane, of the irregular (random) waves on the support structure, and of the turbine re-
sponse. Obtaining realistic response of the turbine depends, among other factors, on
appropriate modeling of the incident wind and waves. The current practice for model-
ing waves on offshore wind turbines is limited to the representation of linear irregular
waves. While such models are appropriate for deep waters, they are not accurate rep-
resentations of waves in shallow waters where offshore wind turbines are commonly
sited. In shallow waters, waves are generally nonlinear in nature. It is, therefore, of
interest to assess the influence of alternative wave models on the behavior of wind
turbines (e.g., on the tower response) as well as on extrapolated long-term turbine
loads. The expectation is that nonlinear (second-order) irregular waves [41, 91, 92]
can better describe waves in shallow waters. In this study, we investigate differences
in turbine response statistics and in long-term load predictions that arise from the
use of alternative wave models.
133
Among several nonlinear irregular wave models available in the literature, the
one that has been recommended by offshore guidelines [17] and that has been increas-
ingly applied to a variety of problems in recent years [17, 27, 44, 102] is the second-
order nonlinear irregular wave model developed by Sharma and Dean [92]. This
second-order nonlinear irregular wave model is based on the solution of the bound-
ary value problem—Laplace equation in velocity potential associated with nonlinear
boundary conditions—using a second-order perturbation expansion of the relevant
variables (velocity potential and sea surface elevation) and a Taylor series expan-
sion of nonlinear boundary conditions about the free surface. Such an approach was
first presented by Longuet-Higgins [57, 58] and Hasselmann [31] to describe nonlinear
irregular waves for infinite water depths. Sharma [91] and Sharma and Dean [92]
and Hudspeth [40, 41] used a similar perturbation approach to derive a nonlinear
irregular wave model for finite water depths. Since its introduction in 1979, this
second-order nonlinear irregular wave model has been studied by several researchers.
Hu and Zhao [39], Langley [53], Longuet-Higgins [59] and Forristall [26], among oth-
ers, investigated the statistical properties (skewness in particular) of this nonlinear
irregular wave model. Jha [44] used this model for simulation of the sea surface
elevation process and compared statistics of simulated waves with those from field
measurements at a site in the North Sea. Forristall [25] compared various stretch-
ing techniques for calculating wave kinematics above the mean sea level. Baldock
et al. [2] compared the deep water approximation of the nonlinear wave model with
experimental results. Due to the maturity of this nonlinear irregular wave model and
because of its theoretical basis (as distinct from some semi-empirical models, such as
the New Wave model [103], the constrained New Wave model [100] and a hybrid wave
model [108]), and also due to its computational efficiency for simulations (as distinct
from realistic though computationally very expensive models based on Boussinesq
134
theory [61, 63, 75]), we will use this second-order nonlinear irregular wave model, as
developed by Sharma and Dean [92], to simulate nonlinear irregular waves in shallow
water depths for offshore wind turbine loads predictions.
The outline of this chapter is as follows: we first describe the theory of the
second-order nonlinear irregular wave model. We limit our focus to unidirectional
(long-crested) waves. We then present a detailed procedure for efficient simulation
of the nonlinear sea surface elevation process and of the water particle kinematics.
We discuss how these theoretical estimates for statistics of this nonlinear wave model
may be computed; these can be used as a check of the accuracy of simulations. Next,
we present two examples to illustrate the nonlinearity in this wave model. We then
compute hydrodynamic loads on a rigid stand-alone monopile, and discuss the me-
chanics of the wave loading and how the nonlinear waves influence hydrodynamic
loads. We incorporate the nonlinear wave model in the turbine response simulation
software, FAST [47], for the aero-servo-hydro-elastic analysis of offshore wind tur-
bines, and we compute loads on the monopile support structure of a NREL 5MW
offshore wind turbine model for several representative environmental states where we
focus on differences in tower bending moment statistics at the mudline due to linear
and nonlinear waves. Finally, we compare long-term load predictions using inverse
FORM with both linear and nonlinear wave models. We discuss convergence crite-
ria in predicting accurate 20-year loads and discuss wind versus wave dominance in
the load prediction. In the following, the term ‘linear waves’ refers to linear irregular
waves and the term ‘nonlinear waves’ refers to second-order nonlinear irregular waves,
unless otherwise indicated.
135
5.2 Second-order Nonlinear Irregular Waves
Linear wave theory for regular or irregular waves involves solution of Laplace’s
equation expressed in terms of a velocity potential and the use of linearized boundary
conditions [90]. For nonlinear waves, the theory involves application of a perturbation
approach to solve Laplace’s equation with nonlinear boundary conditions. Sharma
and Dean [92] and Hudspeth and Chen [41] used such an approach to derive a nonlin-
ear wave theory for finite water depths. We will use the formulation of Sharma and
Dean, which is briefly described below.
The nonlinear sea surface elevation, η(x, t), at a spatial location, x (see Fig. 4.5
for a definition of the coordinate axes), and at a time instant, t, may be expressed as a
sum of the first-order component, η1 (x, t), and the second-order component, η2 (x, t):
The first-order component, η1 (x, t), (earlier discussed in Eq. 4.7 in Chapter 4 where
we denoted it simply as η(x, t), we present it here again only for completeness) is
given as in linear wave theory by
N
X
η1 (x, t) = Am cos(ωm t − km x − φm ) (5.2)
m=1
where ωm refers to the frequency of the mth wave component and φm is the associated
random phase assumed uniformly distributed over [0, 2π]. The variable, km , is the
wave number which is related to the frequency, ωm , through the linear dispersion
relation given as follows:
2
ωm = gkm tanh(km h) (5.3)
136
The amplitudes of the wave components, Am , in Eq. 5.2, are Rayleigh-distributed
random variables whose mean square value is given as follows:
where S(ωm ) refers to the one-sided power spectral density function of the sea surface
elevation process. The integer, m, refers to a frequency index that ranges from 1 to
N , the total number of wave components represented in the simulated wave train.
The second-order component, η2 (x, t), is obtained as a result of the sum- and
difference-frequency interactions as follows:
N X
X N
£ © − +
ª¤
η2 (x, t) = Am An Bmn cos(ψm − ψn ) + Bmn cos(ψm + ψn ) (5.5)
m=1 n=1
ψm = (ωm t − km x − φm ) (5.6)
− +
The second-order transfer functions, Bmn and Bmn , in Eq. 5.5 are obtained from
the solution of Laplace’s equation for the velocity potential with nonlinear boundary
− +
conditions. They (i.e., Bmn and Bmn ) are functions of frequency and wave number
and are independent of the spectrum used. Equations for these transfer functions are
as follows:
· − ¸
− 1 Dmn − (km kn + Rm Rn )
Bmn = √ + (Rm + Rn ) (5.7)
4 Rm Rn
· + ¸
+ 1 Dmn − (km kn − Rm Rn )
Bmn = √ + (Rm + Rn ) (5.8)
4 Rm Rn
137
where
√ √ ©√ 2 2
√ ª
− ( Rm − Rn ) Rn (km − Rm ) − Rm (kn2 − Rn2 )
Dmn = √ √ (5.9)
( Rm − Rn )2 − kmn − tanh(k − h)
mn
√ √
2( Rm − Rn )2 (km kn + Rm Rn )
+ √ √
( Rm − Rn )2 − kmn
− tanh(k − h)
mn
√ √ ©√ 2 2
√ ª
+ ( Rm + Rn ) Rn (km − Rm ) + Rm (kn2 − Rn2 )
Dmn = √ √ (5.10)
( Rm + Rn )2 − kmn+ tanh(k + h)
mn
√ √ 2
2( Rm + Rn ) (km kn − Rm Rn )
+ √ √
( Rm + Rn )2 − kmn
+ tanh(k + h)
mn
− +
The variables, Rm , kmn and kmn , are given as follows:
2
ωm
Rm = (5.11)
g
−
kmn = |km − kn | (5.12)
+
kmn = km + kn (5.13)
The first-order velocity potential, Φ1 (x, z, t), is given by Eq. 4.10. The second-
order velocity potential can be broken into its sum- and difference-frequency interac-
tion terms, such that,
Φ2 (x, z, t) = Φ− +
2 (x, z, t) + Φ2 (x, z, t) (5.15)
N N · ¸
1 XX −
cosh(kmn (h + z)) Dmn−
Φ−
2 (x, z, t) = bm bn − h)
sin(ψm − ψn )
4 m=1 n=1 cosh(kmn (ωm − ωn )
N N · ¸
+ 1 XX +
cosh(kmn (h + z)) Dmn+
Φ2 (x, z, t) = bm bn + h)
sin(ψm + ψn )
4 m=1 n=1 cosh(kmn (ωm + ωn )
138
where
Am g
bm = (5.16)
ωm
The horizontal water particle velocity, u(x, z, t), may be obtained from the
spatial derivative of the velocity potential (as in Eq. 4.12), and the horizontal water
particle acceleration, u̇(x, z, t), may be obtained as the temporal derivative of the
velocity (Eq. 4.13).
The theoretical model for representing the second-order sea surface elevation
and water particle kinematics presented above needs to be implemented numerically
for digital simulation of waves. While this second-order wave model has been discussed
in several published studies [17, 26, 38, 41, 91], a complete recipe for the numerical
simulation of such second-order waves is not yet available in the literature. This is
in contrast to the case for linear irregular waves, for which a recipe for numerical
simulation with all the pertinent details is available in several places (e.g., [3, 93]).
It is important that all the details for efficient numerical simulation of second-order
waves be presented in one place. Therefore, we discuss the numerical simulation of
second-order waves in detail below.
We are interested in simulating the sea surface elevation and particle kinemat-
ics only at the centerline of the monopile, which we arbitrarily choose to be where
x = 0 (see Fig. 4.5), i.e., we simulate η(0, t); henceforth, we will denote η(0, t), η1 (0, t),
and η2 (0, t) simply as η(t), η1 (t) and η2 (t), respectively.
The double sum in Eq. 5.5, to simulate the second-order component of the
139
sea surface elevation process, may be rewritten as
" N N #
XX£ ¤
±
η2 (tp ) = < Xmn exp (−i(m ± n)∆ωp∆t) (5.17)
m=1 n=1
± ±
Xmn = Am An Bmn exp(−i(φm ± φn )) (5.18)
Performing this sum (as given by Eq. 5.17) is computationally inefficient even
for moderate values of N , as a total of N 2 terms are to be summed. Therefore, we seek
an implementation based on the IFFT (inverse Fast Fourier Transform) technique.
We start by rewriting the double summation in Eq. 5.17 as an equivalent single
summation as follows:
" N #
X
η2 (tp ) = < Yj± exp(−i(j∆ω)(p∆t)) (5.19)
j=1
Note that the term within the summation in Eq. 5.20 is the definition of the inverse
discrete Fourier transform of Yj± , which can be efficiently computed using inverse
Fast Fourier transform (IFFT) techniques [76], such that
£ ¡ ¢¤
η2 (tp ) = < IFFT Yj± . (5.21)
The Fourier coefficients Yj± are obtained by “equating” Eq. 5.20 to Eq. 5.17. This
involves conversion of a two-dimensional sum over N 2 terms to a one-dimensional
140
sum over N terms by collecting all (m, n) pairs that yield a sum, or difference,
equal to j, for sum- and difference-frequency interactions, respectively. For exam-
ple, (m, n) = (1, 3), (2, 2) and (3, 1) would all yield a sum for j = 4; so, these three
+
Xmn contributions must be collected to form Y4+ . Note that both the time index, p,
and the frequency indices, j, still range from 1 to N .
Recall that, for simulation of first-order waves, the sampling theorem (Eq. 4.15)
requires that the sampling frequency, ωs , must be at least twice the cut-off frequency,
ωcut ; if ωcut = M ∆ω, then, we must satisfy N ≥ 2M and S(m∆ω) = 0 for M + 1 <
m < N . Therefore, each term in the summation in Eq. 5.17 is non-zero only over the
first M frequencies. As a result, the largest frequency at which the terms in Eqs. 5.17
or 5.20 are non-zero is associated with the frequency index, 2M (for sum-frequency
interactions). Therefore, to satisfy the sampling theorem for second-order simulation,
we must ensure (by simply defining the spectrum in an appropriate way) that
N ≥ 4M (5.22)
This simply means that only the first 2M terms, out of N in Eq. 5.20 are non-zero
and need to be evaluated. Details related to collecting the sum- and the difference-
frequency interaction terms to form Yj± required in Eq. 5.20 are presented next.
Sum-Frequency Interactions
141
number of (m, n) pairs to be summed for a single value of j, then it may be seen that:
0 j<2
j−1 2≤j≤M
r(j) = (5.24)
2M − (j − 1) (M + 1) ≤ j ≤ 2M
0 j > 2M
XX
Yj+ = +
Xmn (5.25)
| {z }
m+n=j
with Y1+ = 0.
Difference-Frequency Interactions
(Eq. 5.20) is the same as the summation-based simulation (Eq. 5.17). Following from
the symmetry of the cosine function, we need
½
(φm − φn ), m>n
φ−
mn = (5.26)
−(φm − φn ), m ≤ n
Again, more than one (m, n) pair may result in a difference, j = |m − n|. If r(j) is
the number of (m, n) pairs to be collected for a single value of j, then it may be seen
that:
½
2(M − j) 0 ≤ j ≤ (M − 1)
r(j) = (5.27)
0 j≥M
142
The Fourier coefficients Yj− for 1 ≤ j ≤ (M − 1) are obtained as
XX
Yj− = −
Xmn (5.28)
| {z }
|m−n|=j
Using the procedure described above, we can obtain the coefficients, Yj± , for
Eq. 5.20 from Eqs. 5.25 and 5.28, and simulate the time series of the second-order sea
surface elevation, η2 (tp ), using Eq. 5.21.
± ±
Umn = Zmn exp(−i(φm ± φn )sgn(m ± n)) (5.31)
± ±
U̇mn = −(ωm ± ωn )Umn (5.32)
± ±
± 1 cosh(kmn (h + z)) Dmn ±
Zmn = bm bn ±
kmn (5.33)
4 cosh(kmn h) (ωm ± ωn )
143
The simulation of particle velocity and acceleration according to Eqs. 5.29
and 5.30, respectively, requires a double summation. We can reduce these double
summations to single summations just as we did for sea surface elevation, and can
readily simulate particle velocity and acceleration according to the following equa-
tions:
£ ¡ ¢¤
u2 (z, tp ) = < IFFT Wj± (5.35)
h ³ ´i
u̇2 (z, tp ) = = IFFT Ẇj± (5.36)
The Fourier coefficients, Wj± , in Eq. 5.35 and Ẇj± in Eq. 5.36, are assembled
± ±
from the coefficients, Umn and U̇mn , respectively, using exactly the same approach as
for the coefficients, Yj± , in Eq. 5.20 that were in turn assembled from the coefficients,
±
Xmn , in Eq. 5.17.
The area under the power spectral density function of the sea surface elevation
process is equal to the variance of a simulated linear sea surface elevation ensemble,
a fact that can be used to evaluate the accuracy of the simulations. Likewise, one
can compute the skewness from wave spectra that should match that from a sim-
ulated nonlinear sea surface elevation ensemble. In fact, nonlinearity described by
the second-order irregular wave model is directly related to the process skewness.
Langley [53] developed a procedure, based on Volterra series models, to compute the
statistical moments of a nonlinear sea surface elevation process, which we describe in
the following.
We first define matrices, Q and P , whose (m, n) elements are qmn and pmn ,
144
and a vector, s, whose elements, sm , are such that
− +
p
qmn = sm sn Bmn , pmn = sm sn Bmn , sm = S(ωm )dω (5.37)
− +
where S(ωm ) is the wave spectrum, and Bmn and Bmn are the second-order transfer
functions given by Eqs. 5.7 and 5.8. We then form matrices, T1 and T2 , such that
T1 = Q + P, T2 = Q − P (5.38)
Langley [53] derived the following expressions for the first four cumulants, cn ,
in terms of β and λ:
2N
X
c1 = λj (5.42)
j=1
2N
X
c2 = 2λ2j + βj2 (5.43)
j=1
2N
X
c3 = 8λ3j + 6βj2 λj (5.44)
j=1
2N
X
c4 = 48(λ4j + βj2 λ2j ) (5.45)
j=1
145
Note that for a zero-mean process, the first three cumulants, c1 , c2 and c3 , are the
same as the first three moments, m1 , m2 and m3 . The fourth moment, m4 , is related
to the fourth cumulant, c4 , as follows:
m4 = c4 + 3m22 (5.46)
Some observations regarding the statistical model presented above are in order.
Note that due to the definition of matrices, T1 and T2 (Eq. 5.38), the trace of the
eigenmatrix, Λ1 , is equal to the negative of the trace of the eigenmatrix, Λ2 . As a
result, the first cumulant, c1 (Eq. 5.42), or the first moment, is identically zero, which
implies that the nonlinear wave elevation process has zero mean. The second term
in Eq. 5.43 is the area under the power spectrum, σ 2 , which is the contribution of
the first-order process, η1 , to the mean square value. The first term in Eq. 5.43 is
the contribution of the second-order process, η2 , to the mean square value and this is
comparatively smaller generally. To compute the variance of the nonlinear sea surface
elevation process, contributions from both the first- and second-order components
should be included. The area under the power spectrum should be compared to the
variance, σ 2 , of the first-order process. The first term in Eq. 5.44 is equal to E[η23 ]
and this is usually negligible (E[·] denotes the expected value operator). The second
term corresponds to E[3η12 η2 ], which dominates the overall skewness of the process.
3/2
The skewness of the nonlinear process, η, is m3 /m2 and not m3 /σ 3/2 . If the latter
formulation is used, a larger skewness would result since m22 > σ 2 . It is important to
note that the theoretical statistics computed using the formulation presented above
would match statistics computed from an ensemble of several time series realizations,
and not from a single time series. Even for linear wave simulations, only the average
of variance computed from several simulations will exactly match the area under the
target wave spectrum. Also, Hu and Zhao [39] pointed out that while the skewness
146
averaged over multiple simulations may match the theoretical value, the coefficient
of variation on this skewness estimate may be quite large. To determine how many
simulations are required to obtain a stable estimate of the average skewness, it might
be necessary to carry out a convergence study.
where Hs is the significant wave height and Lz is the wavelength (obtained using the
linear dispersion relation) corresponding to the zero-crossing frequency, ωz , equal to
p
γ2 /γ0 ; here, γ0 and γ2 , are the zeroth and second spectral moments, respectively.
The nth spectral moment, γn is given by
Z
γn = ω n S(ω)dω (5.48)
Based on results from numerical simulations, Hu and Zhao [39] suggested that the
proposed second-order wave model is valid as long as the wave steepness is smaller
than approximately 0.08. We will use this empirical limit for our simulations.
147
η1 η+2 η−2 η
2
Figure 5.1: Simulation of sea surface elevation using a delta function power spectrum
showing differences between linear and nonlinear waves.
The simplest example that clearly illustrates the difference between linear and
nonlinear waves is a wave simulated using a delta function spectrum. Consider a delta
function spectrum which is one with energy concentrated at a single frequency, ω0
(taken to be 0.5 Hz here). Simulation of a linear irregular wave for this spectrum
effectively results in a regular sinusoidal wave with a period of 2 seconds. Figure 5.1
shows that the sum-frequency component of the second-order wave, η2+ , is a wave
with frequency, 2ω0 . This component, when added to the first-order wave, results in
systematically higher crests and shallower troughs; this is what gives rise to asymme-
try in the waves. For irregular waves, such asymmetry results in a non-zero skewness
of the sea surface elevation process. The difference-frequency component, η2− , in this
case, corresponds to a zero-frequency signal that has zero amplitude; for a spectrum
with energy present at more than one frequency, the difference-frequency component
148
η1
−5
−10
Figure 5.2: Time series of the sea surface elevation process simulated using a JON-
SWAP spectrum with Hs = 12 m and Tp = 14 sec.
typically results in a small set-down effect. For structures with low natural frequen-
cies, the difference-frequency effect may become important.
For real seas, the sea surface elevation spectrum has energy distributed over
a range of frequencies and a parametric spectral representation such as the JON-
SWAP spectrum [17] is often used. We consider now a JONSWAP spectrum with
a significant wave height of 12 m and a peak spectral period of 14 sec. The water
depth is assumed to be 100 m. This example was also used by Sweetman and Win-
terstein [97] to demonstrate the capabilities of their computer program, WaveMaker,
for simulating nonlinear waves (they did not, however, discuss wave kinematics and
forces). Figure 5.2 shows a segment of 200-second duration for a single realization
of the simulated linear and nonlinear sea surface elevation processes. Clearly, larger
crests are observed for the nonlinear waves as one would expect. Table 5.1 shows that
the variance and skewness of the sea surface elevation process, averaged over 1,000
149
Table 5.1: Statistical moments of the linear and nonlinear
sea surface elevation processes simulated for a JONSWAP
spectrum with Hs = 12 m and Tp = 14 sec.
Variance (m2 ) Skewness (-)
Wave model → linear nonlinear linear nonlinear
Target1 8.84 8.94 0 0.17
Average2 9.08 9.17 0 0.16
COV 2 0.14 0.15 - 0.42
1
Computed from the spectrum (see Section 5.2.3).
2
Statistics based on 1,000 simulations.
simulations, match well with the target values computed for this spectrum using the
procedure presented in Section 5.2.3. The variance for the nonlinear waves is larger
than that for the linear waves, and the difference is the contribution of the second-
order waves to the variance of the nonlinear wave process (see Eq. 5.43). The skewness
of the nonlinear waves is estimated to be 0.16, while that of the linear waves, which
describe a Gaussian process, is zero. Note that the non-zero skewness of steep waves,
which is often observed in field measurements, cannot be captured by the linear wave
theory and, thus, the use of a nonlinear wave theory is desirable [44].
150
0.5
Linear waves
0.3
Nonlinear waves
Probability of Exceedance
Gaussian distribution
0.1
0.01
0.001
0.0001
1e−005
1e−006
1e−007
0 4 8 12 16
Sea surface elevation (m)
due to the positive skewness (and also due to the larger variance), the sea surface
elevation for low (rare) exceedance probabilities is larger for nonlinear waves than it
is for linear waves. In fact, the largest observed sea surface elevation for linear waves
is 13.8 m while that for nonlinear waves is 15.6 m — about 14% larger. The average,
over all simulations, of the maximum sea surface elevation from each simulation, is
10.3 m for nonlinear waves, which is about 10% larger than that (9.3 m) for linear
waves.
An interesting point worth noting in Table 5.1 is the rather large coefficient of
variation (COV) of 41% on the estimate of the skewness value. Such large variation
is expected as is discussed by Hu and Zhao [39]. In fact, this is the reason why a large
ensemble of random realizations is required if higher moments computed from time
series are to match target values. To illustrate that a larger number of simulations
are required to obtain a stable estimate of the process skewness than that required to
151
obtain a stable estimate of the process variance, we carry out a convergence check on
an error estimate. If the average of calculated values of a statistic, θ, over N samples,
θ̂N , is an estimate of θ, and we assume that the average over a large (1,000 here)
sample size, θ̂∞ , is a stable estimate, then we can define a measure of error, ² (in
percent), as follows.
|θ̂∞ − θ̂N |
²= × 100 (5.49)
θ̂∞
We can now find the number of simulations required to obtain a desired error.
It was found, in our studies, that 50 simulations are required to obtain a stable
estimate of the variance, while 410 simulations are required for a stable estimate of
the skewness if the error, ², is to be at most 3%. For a stricter error bound of 1%,
260 and 740 simulations are required for variance and skewness, respectively. This
clearly shows that for any specified error criterion, a larger number of simulations is
required to obtain a stable estimate for skewness than is required for variance of the
sea surface elevation process.
152
waves, while using the Wheeler stretching technique to compute kinematics up to the
instantaneous sea surface. We use Morison’s equation [3] to compute hydrodynamic
loads per unit length, f . Thus, we have
1 πD2
f = fD + fM = CD ρDu|u| + CM ρu̇ (5.50)
2 4
where fD and fM are drag and inertia forces, respectively, CD is the drag coefficient
taken as 1.0, CM is the inertia coefficient taken as 2.0. These are properties assumed
for the NREL 5MW turbine model [48]. In Eq. 5.50, u is the horizontal water particle
velocity, u̇ is the horizontal water particle acceleration, ρ is the density of water, and
D is the monopile diameter. Note that Eq. 5.50, unlike Eq. 4.30, does not include the
effect of the relative motion between water particles and the structure, because we are
only carrying out a static analysis here. The sea state, for the following discussions,
is the one that governed long-term loads for the turbine (discussed in Table 4.3) when
linear waves were modeled, with a significant wave height of 5.5 m and spectral peak
period of 11.2 sec.
153
5.3.1 Comparison of Loads due to Linear and Nonlinear Waves
Table 5.2 shows that the maximum quasi-static load due to nonlinear waves
on the monopile is 36.5 MN-m, which is about 29% larger than the corresponding
maximum load of 28.3 MN-m due to linear waves; thus, nonlinear waves are clearly
important. Moreover, the standard deviation, skewness, kurtosis, and peak factor are
all larger with the nonlinear waves; differences in these statistics indicate a strongly
non-Gaussian character introduced by the nonlinear waves to the load process. Large
kurtosis values (greater than 3), for instance, imply fatter tails than for a Gaussian
distribution; this means that larger loads might be more likely to occur. The larger
peak factors on ten-minute maximum loads seen with nonlinear waves are evidence of
this and are indicative of maximum values systematically at higher levels above the
mean value. What all of this suggests is that the ratio of extrapolated loads based
on the use of nonlinear waves to that based on linear waves could be larger than the
corresponding ratio of 29% for maximum loads before extrapolation.
60
1000
40
500
20
0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Frequency (Hz) Frequency (Hz)
(a) (b)
Figure 5.4: Power spectrum of (a) sea surface elevation and (b) bending moment
at the base of the rigid stand-alone monopile, when linear and nonlinear waves are
simulated using a JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec.
154
Figure 5.4(a) shows the power spectral density (PSD) for the sea surface ele-
vation process. The PSDs for both linear and nonlinear waves have a main peak at
about 0.09 Hz corresponding to the spectral peak period of the waves which is 11.2
sec. For nonlinear waves, there is a small second peak at about 0.18 Hz, which is about
twice the spectral peak frequency of the waves; it arises due to the sum-frequency
interactions that occur in nonlinear waves. The PSD for nonlinear waves also shows
a small peak at almost zero frequency, which is due to the difference-frequency inter-
actions in nonlinear waves. The power spectra of the bending moment at the base
of monopile, shown in the Fig. 5.4(b), depicts similar behavior—a secondary peak
at twice the spectral peak period for nonlinear waves. For a dynamic analysis, such
secondary peaks may become important if they are close to the natural period of the
structure.
155
Table 5.3: Comparison of hydrodynamic loads on a stand-alone
monopile analyzed quasi-statically and water particle kinematics,
averaged over fifty simulations, computed with linear and nonlin-
ear irregular waves using a JONSWAP spectrum with Hs = 5.5 m
and Tp = 11.2 sec.
Bending moment at base (MN-m) umax u̇max
Drag Inertia Total (m/s) (m/s2 )
Linear 8.6 26.8 28.3 3.8 3.0
Nonlinear 13.1 32.3 36.5 4.8 3.5
Ratio 1.52 1.21 1.29 1.25 1.20
Note: umax and u̇max are maximum water particle velocity and
acceleration, respectively, computed at the free surface.
Another observation from Table 5.3 is that the total hydrodynamic loads on
this monopile are dominated by inertia forces, both for linear as well as nonlinear
waves. In the linear case, the ratio of the maximum tower base moment due to inertia
forces to that due to drag forces is more than three; in the nonlinear case, this ratio
is about 2.5 (again, this is because of the greater influence of nonlinear waves on drag
forces compared to inertia forces). For linear waves, whether a structure is drag- or
inertia-dominated may be evaluated without carrying out time-domain simulations.
To study the relative importance of drag versus inertia, one can examine the ratio of
the standard deviations of the drag load to the inertia load (since both are random
processes); this has been studied by others [18, 62] and the ratio can conveniently be
written as follows:
r
σ[fD ] 3
= κ (5.51)
σ[fM ] 4
CD Hs cosh k(z + h)
κ ≈ (5.52)
πCM D sinh(kh)
where κ is analogous to the Keulegan-Carpenter number [18] for regular linear waves;
also, h is the water depth and z refers to a vertical coordinate measured positive
156
upwards and with origin at the mean sea level. Once a wave spectrum is given, one
can evaluate Eq. 5.51 for a wave number corresponding, say, to the spectral peak
period where wave energy is high. Since drag forces are maximum at crests, one can
then obtain the maximum value of the ratio of the drag to inertia forces by assuming
z = Hmax /2, where Hmax is the maximum instantaneous wave height which can be
based on estimates of extremes of a Gaussian process representing the sea surface
elevation. A very conservative estimate [62] of Hmax is 2Hs which we use here. With
these considerations and by using dimensions and other parameters specified for the
monopile, the ratio in Eq. 5.51 is found to be approximately 0.2, which confirms that
the monopile under consideration is dominated by inertia forces. We again mention
that loads on this inertia-dominated monopile increased by about 29% when nonlinear
waves are used. Nonlinear waves would have even more pronounced effect on loads
for a drag-dominated structure, such as a jacket structure with slender members.
157
Table 5.4: Comparison of statistics, averaged over fifty simulations, of bending mo-
ment at the base of the rigid stand-alone monopile due to drag and inertia forces,
when linear and nonlinear waves are simulated using a JONSWAP spectrum with Hs
= 5.5 m and Tp = 11.2 sec.
sian factor; therefore, drag forces are non-Gaussian. Both the drag and inertia forces
become more non-Gaussian for nonlinear waves, because nonlinear waves themselves
are non-Gaussian. We next discuss theoretical reasons for the non-Gaussian character
of drag forces for loading due to linear waves.
f = fD + fM = KD u|u| + KM u̇ (5.53)
where
1 πD2
KD = CD ρD; KM = CM ρ (5.54)
2 4
158
where fD and fM are drag and inertia forces, respectively; all the symbols appearing
in Eqs. 5.53 and 5.54 were explained in Eq. 5.50.
A linear irregular sea surface elevation process is Gaussian. The water particle
velocity and acceleration processes, at any depth, simulated using their respective
spectra, Su (ω) and Su̇ (ω) (Eqs. 4.28 and 4.29), are also Gaussian. In fact, it has
been shown that any process simulated using a power spectrum and random phases
is Gaussian [4].
where σu̇ is the standard deviation of the horizontal water particle acceleration ran-
dom process. The probability distribution function (PDF) of Y , fY (y), has been
shown to be given as follows:
r µ 2¶Z ∞ µ ¶
α y − 12 t2
fY (y) = exp − t exp −αt − cosh(yt)dt (5.56)
2π 2 2 0 2
159
The variable, σu , is the standard deviation of the horizontal water particle velocity
random process.
All moments of the random variable, Y , can be computed using its probability
distribution. Since fY (y) is an even function, all odd moments of Y are zero. The
second and fourth moments are given by:
3
E[Y 2 ] = +1 (5.58)
4α2
105 9
E[Y 4 ] = 4
+ 2 +3 (5.59)
16α 2α
E[Y 4 ]
γ= (5.60)
E[Y 2 ]2
160
Table 5.5: Comparison of statistics of Morison forces (per unit length) at
the mean sea level, averaged over fifty simulations, with statistics obtained
from the theoretical model, when linear irregular waves are simulated with
a JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec.
Standard deviation (kN/m) Kurtosis
Hydrodynamics → drag inertia total drag inertia total
Simulations 7.0 49.1 49.6 10.7 3.0 3.0
Theoretical 7.4 49.2 49.7 9.9 3.0 3.0
character of drag forces. Since the bending moment at the base of the monopile is
obtained by integrating the Morison forces at all nodes, it (the bending moment) is
also expected to be non-Gaussian with large kurtosis values, as has been shown by
Najafian and Burrows [72]; this is the reason we observed such large kurtosis values
for the bending moment at the base of the monopile in Table 5.4. Note, however,
that the Morison forces in Table 5.5 are based on unstretched wave kinematics, while
the bending moments in Table 5.4 are computed using stretched kinematics, which
introduces another source of nonlinearity. While the above discussion illustrates the
non-Gaussian character of Morison forces due to linear irregular waves, the non-
Gaussian character of forces only increases (as shown in Table 5.4) in presence of
nonlinear irregular waves, which are inherently non-Gaussian.
For a given significant wave height, as the spectral peak period is increased, the
corresponding spectral peak frequency and the wave number decrease. Consequently,
the wavelength increases and the wave steepness decreases. Wave steepness has a
direct influence on the particle kinematics and therefore on Morison forces. It is of
interest to investigate how hydrodynamic loads computed from linear and nonlinear
161
Table 5.6: Values of steepness, variance and skewness of the
nonlinear sea surface elevation process for various choices
of spectral peak period, Tp , for a significant wave height,
Hs , of 5.5 m, when waves are simulated using a JONSWAP
spectrum.
Sea Surface Elevation
√
Tp / Hs Tp (sec) steepness variance (m2 ) skewness
5.0 11.7 0.050 2.00 0.11
4.4 10.3 0.061 1.97 0.15
4.0 9.4 0.071 1.95 0.18
3.6 8.4 0.084 1.93 0.21
In the simulations discussed thus far (Tables 5.2-5.4), for the governing sea
state (Table 4.3) with a significant wave height, Hs , of 5.5 m, we have used only one
value of the spectral peak period, Tp , of 11.2 sec, which corresponds to a steepness
value of 0.049. This choice was based on an empirical correlation between Hs and
Tp for the chosen site (NDBC buoy 44007). A more general range of spectral peak
periods, for the JONSWAP spectrum, is suggested in the DNV guidelines [17], which
is given as follows:
Tp
3.6 < √ < 5.0 (5.61)
Hs
162
Bending moment at the base of the rigid stand−alone monopile (MN−m)
40
35
30 Linear waves
Nonlinear waves
25
0.05 0.06 0.07 0.08 0.09
Ratio of bending moment from nonlinear waves to that from linear waves
1.4
1.3
1.2
1.1
1
0.05 0.06 0.07 0.08 0.09
Wave steepness
Figure 5.5: Variation of the ten-minute maximum bending moment at the base of
a rigid stand-alone monopile, averaged over fifty simulations, with wave steepness,
when linear and nonlinear waves are simulated using a JONSWAP spectrum with Hs
= 5.5 m while Tp values are chosen according to the wave steepness.
peak period. Table 5.6 shows that the wave steepness increases as the wave period
is decreased. Since the waves are more asymmetric for larger steepness, the skewness
is also seen to increase. The variance of the nonlinear sea surface elevation process
is rather insensitive to wave period; the variance of the linear sea surface elevation
process, which is the area under the sea surface elevation spectrum, is completely
independent of wave period as it depends only on significant wave height.
Figure 5.5 shows the variation of bending moment, at the base of a rigid stand-
alone monopile, with wave steepness. It is found that loads increase with steepness
both for linear and nonlinear waves. However, as the steepness is increased, the loads
increase by relatively larger amounts for linear waves; consequently, the ratio of loads
from nonlinear waves to loads from linear waves decreases with increasing steepness.
163
Water particle acceleration (m/s2)
5
3 Linear waves
Nonlinear waves
2
0.05 0.06 0.07 0.08 0.09
Ratio of particle acceleration from nonlinear waves to that from linear waves
1.3
1.2
1.1
1
0.05 0.06 0.07 0.08 0.09
Wave steepness
Figure 5.6: Variation of the ten-minute maximum horizontal water particle accelera-
tion at the mean sea level, averaged over fifty simulations, with wave steepness, when
linear and nonlinear waves are simulated using a JONSWAP spectrum with Hs = 5.5
m while Tp values are chosen according to the wave steepness.
Since inertia forces dominate Morison forces for this monopile, and inertia
forces at a node are directly proportional to the water particle acceleration at that
node, we study the variation of the horizontal water particle acceleration at the mean
sea level (as acceleration decreases with water depth); these accelerations are studied
versus wave steepness in Fig. 5.6. As was the case with monopile loads, the water
particle acceleration increases, both with linear and nonlinear waves; also, the ratio
of the particle acceleration for nonlinear waves to that for linear waves decreases with
increasing steepness.
To seek a theoretical explanation for the observed behavior of forces and ac-
celerations with variation in wave steepness, we focus on waves for a delta spectrum,
which effectively represents regular waves. Our aim is to derive analytical expressions
164
for amplitudes of the water particle acceleration for linear and nonlinear waves, and
to study how they vary with the wave steepness. We assume that the delta spectrum
has all its energy, S, represented at a single frequency, ω, or equivalently, at the
period, T , where T = 2π/ω.
Let us denote the amplitude of the water particle acceleration at the mean
sea level for linear (first-order) waves as a1 , and that for second-order waves as a2 .
The amplitude of the water particle acceleration for nonlinear waves is a = a1 + a2 .
We wish to investigate accelerations for steepness values, s1 and s2 , respectively, such
that s2 > s1 ; the corresponding frequencies are ω1 and ω2 such that ω1 > ω2 . We
are interested in determining, if the behavior observed in Fig. 5.6 is correct, namely
if the following holds:
a1,s1 + a2,s1 a1,s2 + a2,s2
> (5.62)
a1,s1 a1,s2
where a1,s1 and a1,s2 are amplitudes of the first-order water particle accelerations for
steepness values of s1 and s2 , respectively. Likewise, a2,s1 and a2,s2 are amplitudes
of the second-order water particle accelerations for steepness values of s1 and s2 ,
respectively. After some rearrangement, the inequality above may be equivalently
written as follows:
a1,s2 a2,s2
> (5.63)
a1,s1 a2,s1
The amplitude of the water particle acceleration at the mean sea level for linear waves,
a1 , can be derived from Eq. 4.23, and is given as follows:
a1 = Agk (5.64)
p
where A = S(ω)dω is the wave amplitude, which is kept constant here, g is the
acceleration due to gravity, and k is the wave number corresponding to the frequency,
165
ω. Therefore, the ratio of the amplitudes of the linear water particle accelerations for
two steepness values is simply the ratio of the corresponding wave numbers.
a1,s2 k2
r1 = = (5.65)
a1,s1 k1
The amplitude of the water particle acceleration at the mean sea level for the
second-order waves, a2 , which can be derived from Eq. 5.30, is given as follows.
A2 g 2 k +
a2 = D (5.66)
2 ω2
where D+ is the second-order transfer function (Eq. 5.10) for sum-frequency inter-
actions at (ω + ω), i.e., at 2ω. The ratio of amplitudes of the second-order particle
accelerations is, therefore, given as:
µ ¶2
a2,s2 k2 ω1 D2+
r2 = = (5.67)
a2,s1 k1 ω2 D1+
According to Eq. 5.63, we are interested in showing that the ratio, r1 /r2 , is
greater than unity. Using Eqs. 5.65 and 5.67, we can write this ratio as follows:
µ ¶2
r1 ω2 D1+
= (5.68)
r2 ω1 D2+
Here, ω2 /ω1 > 1 since ω2 > ω1 . The variation of the second-order transfer function,
D+ , for sum-frequency interactions, (ω + ω), for a frequency, ω, is shown in Fig. 5.7,
from which it is clear that the ratio, αD , equal to D1+ /D2+ , is greater than unity for
ω2 > ω1 . Therefore the ratio, r1 /r2 , is also greater than unity. Note that the ratio of
frequencies is the same as the ratio of steepness values; hence, we can rewrite Eq. 5.68
in terms of the steepness ratio, s2 /s1 , where s2 > s1 .
µ ¶2
r1 s2
= αD > 1 (5.69)
r2 s1
166
−3
x 10
9
8.5
D+ (m−2)
8
7.5
This proves the hypothesis presented in Eq. 5.62. That is, we have shown
that the ratio of the water particle accelerations for nonlinear waves to that for linear
waves decreases as wave steepness (or spectral peak frequency) is increased. While we
have analytically proven this behavior only for the case of a delta function spectrum
(i.e., for regular waves), the rationale for making similar conclusions for irregular
waves is generally based on related arguments as were used for regular waves, at least
qualitatively [18]. This, then, can help explain the results from simulations based on
a JONSWAP spectrum that were presented in Figs. 5.5 and 5.6.
In summary, we have shown that (1) larger wave steepness values result in
larger hydrodynamic loads for both linear and nonlinear waves; and (2) the ratio of
forces from nonlinear waves to those from linear waves decreases as wave steepness is
increased. In other words, for a given significant wave height, loads are large when
steepness is high but the difference in loads from linear and nonlinear wave models is
smallest for such large wave steepnesses.
167
Table 5.7: Ratio of the maximum quasi-static base
moment on a stand-alone monopile computed using
nonlinear waves to that using linear waves (based
on 50 ten-minute simulations).
Site Load Ratio
Hs (m) h (m) Tp (sec) D=6m D=3m
4 10 9.1 1.64 1.79
4 20 8.0 1.01 1.02
8 20 12.9 1.67 1.87
8 30 11.8 1.15 1.23
8 40 11.3 1.02 1.07
12 30 12.7 1.40 1.57
12 40 14.7 1.26 1.33
16 40 14.7 1.48 1.52
Hs : Significant wave height, h: water depth, Tp :
Spectral peak period, D: Monopile diameter.
Thus far, we have studied in detail the response of a 6-meter monopile sub-
jected to a single sea state. It is of interest to study the influence of wave modeling
assumptions for other simulations as well. A few parameters that are of interest to
study are monopile diameter, water depth, and significant wave height. Since we are
interested in loads on offshore wind turbines, we chose representative values of these
parameters based on a quick survey of existing offshore wind farms from which it was
found that monopile diameters of the order of 3 to 6 meters and water depths ranging
from 5 to 30 meters were common. For each water depth, we determined a maximum
possible significant wave height, Hs , before waves would likely break; only Hs values
smaller than this maximum level were considered. The associated spectral peak pe-
riod, Tp , required for the JONSWAP spectrum, is selected to be the one for which the
wave steepness, s = Hs /Lz , is about 0.06, which is a rather severe value of steepness.
168
Table 5.7 shows the ratio of the maximum quasi-static base moment (averaged over
50 ten-minute simulations) for a stand-alone monopile analyzed quasi-statically using
nonlinear waves to that using linear waves.
Table 5.7 shows that monopile loads due to nonlinear waves can, in some cases,
be more than 80% larger than loads due to linear waves. Clearly, for such sites, linear
waves would give unconservative estimates of loads. Several other observations can
be made from this table. For any site, the increased level of load due to the use of
nonlinear waves is larger for a smaller diameter monopile. This is because smaller
members are relatively more drag-dominated, where nonlinear waves are also more
important. For a given wave height, the influence of nonlinear waves decreases with
an increase in water depth. Therefore, modeling of nonlinear waves is important
mainly at shallow-water sites where most offshore wind farms are typically located.
169
Table 5.8: Ten-minute maximum fore-aft tower bending moment (in MN-m) at mud-
line, averaged over 20 simulations, computed with linear and nonlinear irregular waves
using a JONSWAP spectrum with Hs = 5.5 m and Tp = 11.2 sec, and V = 16 m/s.
Ten-minute max fore-aft tower
bending moment at mudline (MN-m)
without wind with wind
Linear 25.7 89.4
Nonlinear 29.9 90.9
Ratio 1.17 1.02
2.0, respectively, for this NREL 5MW offshore wind turbine; the turbine model was
described earlier in Section 4.3.1.
When a linear theory for waves was used for this turbine and site, a mean wind
speed, V , of 16 m/s and a significant wave height, Hs , of 5.5 m were found to represent
the environmental state that governed the long-term tower bending moment at the
mudline for a return period of 20 years, as has been discussed earlier (see Table 4.3).
We now investigate how much tower loads change when nonlinear irregular waves are
employed, for the same environmental conditions.
Table 5.8 shows ten-minute maxima of the fore-aft tower bending moment at
mudline computed using linear and nonlinear waves. When wind is not considered,
loads computed based on modeling the nonlinear waves are about 17% larger than
those due to linear waves, which indicates the importance of nonlinear waves. How-
ever, loads resulting from linear and nonlinear waves are only slightly different when
wind is included and the turbine is in operation. This is because for this (V, Hs ) com-
bination, loads due to waves alone make up only a small fraction (about one-fourth)
170
of the loads due to both wind and waves. Based on this, we might conclude that
20-year loads may not differ greatly whether linear or nonlinear waves are considered.
However, the governing environmental state for long-term loads (computed using in-
verse FORM, for example) may itself change when the wave model is changed. It is
therefore important that we study the influence of wave model choice on loads for
other possible environmental states as well.
It is also interesting to note that for linear waves, the maximum load on the
monopile with turbine, which is about 26 MN-m (Table 5.8), is very close to the load
of about 28 MN-m obtained for the stand-alone monopile analyzed quasi-statically
(Table 5.2). However, the two types of analyses are different. While qualitative
conclusions about differences for linear and nonlinear waves would be similar for the
quasi-static analysis of the stand-alone monopile and the dynamic analysis of the
monopile with turbine, quantitative results are not expected to be same. As an
example, on the stand-alone monopile, the loads from nonlinear waves were about
29% larger than those from linear waves. For the monopile with turbine analyzed
here, the loads are still larger for nonlinear waves, but by only about 17%.
171
Table 5.9: Comparison of statistics of the sea surface elevation process simulated
using linear and nonlinear irregular wave models for a JONSWAP spectrum with Hs
= 7.5 m and Tp = 12.3 sec.
to describe its variability (see Section 4.2.2 for more details). For the selected signifi-
cant wave height, we assume a peak spectral period, Tp , of 12.3 sec, which corresponds
to a wave steepness, s = Hs /Lz , of 0.06. Here, Lz is the wavelength corresponding to
the zero-crossing period, Tz , based on the linear dispersion relation; also Tz is related
to the peak spectral period, Tp , for a JONSWAP wave spectrum [17]. Nonlinear waves
that result from use of the second-order model presented are assumed valid up to a
steepness of around 0.08 [39]. Therefore, a steepness of 0.06 is a case of fairly severe
wave nonlinearity.
172
Table 5.10: Comparison of statistics, averaged over 50 simulations, of the fore-aft
tower bending moment at the mudline for linear and nonlinear waves, for a JONSWAP
spectrum with Hs = 7.5 m and Tp = 12.3 sec, and wind speed, V = 16 m/s.
nonlinear waves. Because of the larger skewness and kurtosis and, hence the somewhat
larger peak factor, associated with nonlinear waves, extremes based on the use of a
nonlinear wave theory are larger than those based on linear waves. Since nonlinear
waves tend to have sharper crests, the maximum of the sea surface elevation is larger
(by about 1 m or 18%) than is the case for linear waves. As a result of these larger wave
heights, it is expected that a greater portion of the monopile would get submerged
and, as a result, the monopile would experience greater lateral base forces due to
hydrodynamic loads. The loads are further amplified because particle velocities and
accelerations are also larger for nonlinear waves (see Table 5.3).
Table 5.10 shows statistics of the fore-aft tower bending moment at the mudline
173
(averaged over 50 ten-minute simulations). When wind is not included, the maximum
load (with both drag and inertia forces accounted for) with the nonlinear waves is
about 54 MN-m, which is about 52% larger than that with the linear waves. When
winds are included, the maximum load with nonlinear waves is about 108 MN-m, only
about 11% larger than that with linear waves—still a significant increase. Table 5.10
also shows that the mean value of the tower bending moment is close to zero with
waves alone as input; this is because wave forces, in the absence of currents, are
zero-mean processes, while wind-induced forces cause a non-zero mean that results
from a non-zero mean longitudinal wind speed. To understand the influence of wave
nonlinearity, we focus our attention on loads in the absence of winds, where it is also
useful to study loads due to drag and inertia forces separately. Table 5.10 shows clearly
that this monopile is dominated by inertia loads (see Section 5.3.2 for more details on
why inertia loads dominate), as is typical for such large-diameter monopile cylinders in
shallow waters. Loads due to inertia forces alone increase from about 31 MN-m (with
linear waves) to 44 MN-m (with nonlinear waves) representing an increase of about
40%. On the other hand, loads due to drag forces alone increase from about 14 MN-m
(with linear waves) to about 26 MN-m (with nonlinear waves)—an increase of nearly
90%. Clearly, wave nonlinearity has a greater influence on drag forces. Note that
nonlinearity of the waves increases both particle velocity and acceleration by similar
amounts, as we showed in Table 5.3; however, since drag forces are proportional to
the square of the particle velocity and inertia forces are only linearly proportional to
the particle acceleration (Eq. 5.50), the increase in loads due to wave nonlinearity is
more significant for drag forces.
174
Gaussian is related to the extent by which its skewness deviates from zero and its
kurtosis deviates from 3.0. Using the mean-upcrossing rate of a random process, we
can compute a theoretical peak factor for a specified exposure time assuming the
process were Gaussian [62], and compare that to the actual peak factor computed
directly and empirically from realizations of the process. A process with a positive
skewness and a kurtosis greater than 3.0 (as is the case here) will generally lead to
larger peak factors. This will, in turn, result in larger extremes associated with any
specified rare probability level compared to those predicted for a Gaussian process.
Note that it is such low exceedance probability levels that are of eventual interest
in predicting long-term loads for ultimate limit states. Table 5.10 shows that skew-
ness, kurtosis and peak factor estimates from the simulations are always larger with
nonlinear waves. As a result, extreme loads predicted based on the use of nonlinear
waves will also be larger. Furthermore, deviations from the Gaussian are greater for
drag forces and, therefore, the influence of nonlinear waves on loads would be even
more pronounced for a structure that is dominated by drag forces—e.g., for slender
members such as those in a jacket structure.
175
Linear Nonlinear
Wave Elevation (m)
5
0
−5
250 300 350 400 450
Wind Speed (m/s)
20
10
PSD [(m)2/Hz]
80
1
10
60
0 40
10
20
−1
10 0
0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4
Frequency (Hz) Frequency (Hz)
(b)
Figure 5.8: (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for wind speed (for V = 16 m/s) and of sea surface
elevation (for linear and nonlinear waves simulated using a JONSWAP spectrum with
Hs = 7.5 m and Tp = 12.3 sec).
176
Linear Nonlinear
Fore−aft Tower Base Shear (MN−m)
4
2
0
−2
−4
250 300 350 400 450
FATBS FATBM
15
5000
PSD [(MN−m)2/Hz]
4000
PSD [(MN)2/Hz]
10
3000
5 2000
1000
0 0
0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4
Frequency (Hz) Frequency (Hz)
(b)
Figure 5.9: (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for fore-aft tower base shear (FATBS) and tower
bending moment (FATBM) at the mudline (for Hs = 7.5 m and Tp = 12.3 sec, when
wind is not included).
177
Linear Nonlinear
Fore−aft Tower Base Shear (MN−m)
4
2
0
−2
−4
250 300 350 400 450
FATBS FATBM
15
5000
PSD [(MN−m)2/Hz]
4000
PSD [(MN)2/Hz]
10
3000
5 2000
1000
0 0
0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4
Frequency (Hz) Frequency (Hz)
(b)
Figure 5.10: (a) A 200-second segment of a representative time series and (b) power
spectral density (PSD) functions for fore-aft tower base shear (FATBS) and tower
bending moment (FATBM) at the mudline (for Hs = 7.5 m and Tp = 12.3 sec, when
wind is included, i.e., V = 16 m/s).
178
To further investigate the influence of wave nonlinearity, we study a represen-
tative 200-second segment from a single ten-minute simulated time history. Figure 5.8
shows that crests of the sea surface elevation process are systematically higher for the
nonlinear wave model than for the linear model. Consequently, positive maxima of
the tower base shear and bending moment at the mudline are also larger for nonlinear
waves, when winds are not included as shown in Fig. 5.9. The power spectrum of the
sea surface elevation (Fig. 5.8) for the nonlinear waves case shows a secondary peak
at about 0.16 Hz, which is twice the spectral peak frequency of 0.08 Hz (since Tp =
12.3 sec), and a small peak close to a zero frequency. Such secondary peaks arise due
to sum- and difference-frequency interactions according to the second-order nonlin-
ear wave model (see Eq. 5.5). These secondary peaks also appear in power spectra of
both the mudline tower base shear and the mudline tower bending moment (Fig. 5.9).
The tower bending moment power spectrum has a significant peak at around 0.27 Hz,
which is the natural frequency for the first bending mode of vibration of the tower in
the fore-aft direction. This clearly suggests that tower dynamics are important in the
overall response. The response of the turbine tower when winds are included is shown
in Fig. 5.10. The influence of winds on the mudline tower base bending moment is
more pronounced than on the mudline tower base shear, because of the large lever
arm associated with the moment due to wind forces acting at the level of the rotor
plane. Comparing time series of the tower bending moment when winds are included
(Fig. 5.10) to when winds are not included (Fig. 5.9), we see that winds add a signifi-
cant mean component to the tower bending moment. The turbulent character of the
wind is also evident. More interesting is the observation that the peak at the tower
natural frequency in the power spectrum of the tower bending moment has almost
disappeared when winds are included. Such behavior has also been reported by other
researchers [46, 101]; it is due to aerodynamic damping from wind loads when the
179
turbine is in operation, which we briefly discuss next.
Our interest is in predicting extreme loads associated with rare (low) proba-
bility of exceedance levels. Earlier (in Section 4.6), with linear irregular waves, we
discussed the subject of statistical extrapolation of wind turbine loads in some de-
180
Table 5.11: Largest observed ten-minute maximum fore-aft tower bending moment
at the mudline, from 500 simulations, for the governing environmental state for the
linear wave model (V = 16 m/s, Hs = 5.5 m) and for the nonlinear wave model (V
= 18 m/s, Hs = 7.5 m).
Ten-minute maximum fore-aft tower
bending moment at mudline
V Hs wave model
(m/s) (m) Linear Nonlinear ratio
(MN-m) (MN-m) -
16 5.5 116.5 117.8 1.01
18 7.5 120.9 140.1 1.16
tail, and derived long-term loads for our 5MW offshore wind turbine model for a
return period of 20 years. We are interested in performing similar calculations, using
3D inverse FORM, for turbine loads using nonlinear waves. To this end, we perform
multiple simulations of the turbine response for various possible (V, Hs ) combinations,
followed by estimation of the fractile, p3 (Eq. 4.35), of load extremes from short-term
distributions of load conditional on environment. The maximum value of the p3 frac-
tile over all environmental states gives the desired long-term load. Such calculations
show that a mean wind speed, V , of 18 m/s and a significant wave height, Hs , of 7.5
m govern the long-term 20-year load when a nonlinear wave model is used. This is
different from the governing values for V of 16 m/s and for Hs of 5.5 m when a linear
wave model was used. Before presenting the actual long-term 20-year loads, we first
compare these two environmental states in the following.
181
0 0
10 10
Exceedance Probability
Exceedance Probability
in 10−minutes
in 10−minutes
−2 −2
10 10
−4 −4
10 10
0 0
10 10
Exceedance Probability
Exceedance Probability
in 10−minutes
in 10−minutes
−2 −2
10 10
−4 −4
10 10
simulations. For comparison, we also ran 500 simulations for V = 16 m/s and Hs =
5.5 m, although this environmental state requires a rarer (3.87 × 10−6 ) probability of
exceedance of load given environment. Table 5.11 shows that the observed maximum
loads in 500 simulations are almost the same with linear and nonlinear wave models
182
for V = 16 m/s and Hs = 5.5 m; this was discussed in detail in Section 5.4.1. For V
= 18 m/s and Hs = 7.5 m, observed maximum loads increase significantly, by about
16%, when the wave model is changed from linear to nonlinear. This is in line with
our earlier discussion (in Section 5.4.2) where we saw that for larger significant wave
heights, the nonlinear wave model results in a significant increase in hydrodynamic
loads. Furthermore, for our pitch-controlled wind turbine, mean loads due to wind
only reduce as the wind speed is increased above rated (11.5 m/s here). Due to a
combination of these two factors, the environmental state of V = 18 m/s and Hs =
7.5 m governs long-term loads when a nonlinear wave model is used.
While we have performed 500 ten-minute simulations for these two environ-
mental states, in practice, a much smaller number of simulations is typically carried
out for extrapolation. We earlier proposed a convergence criterion (Eq. 3.10), which
183
Table 5.12: Estimates of the 84th percentile 10-min maximum fore-aft tower bending
moment at the mudline from 10 separate sets of 50 simulations each.
Table 5.13: Estimates of the 90% relative confidence bound (RCB) on the 84th per-
centile 10-min maximum fore-aft tower bending moment at the mudline.
compared confidence bounds (90% there) on an observed rare percentile (the 84th
percentile there) of the short-term load distribution versus an imposed acceptable
184
limit, to establish how many simulations would be adequate in order to have a sta-
ble distribution especially in tails. Table 5.12 shows estimates of the 84th percentile
10-minute maximum load (fore-aft tower bending moment at the mudline) from ten
separate sets of 50 simulations each. The average of the ten estimates, shown in Ta-
ble 5.12, is very close to the estimate of the 84th percentile load from 500 simulations,
which is 100 MN-m (with linear waves) and 100.7 MN-m (with nonlinear waves) for
V = 16 m/s, Hs = 5.5 m; and 96.1 MN-m (with linear waves) and 108.3 MN-m (with
nonlinear waves) for V = 18 m/s, Hs = 7.5 m. This suggests that 50 simulations may
be sufficient to yield a stable estimate of a rare load fractile. It can also be seen that
the 84th percentile load estimate for all the cases is very stable, and it has a range
(as defined in Table 5.12) that is always less than 10%. This range on the 84th per-
centile loads, which is essentially the 82% (10/11 - 1/11 = 9/11 = 0.82) confidence
interval relative to the median estimate, is larger for nonlinear waves for both the
environmental states. This is expected since nonlinear waves are more non-Gaussian
and highly variable, and hence result in larger variability in turbine loads as well.
The 90% relative confidence bounds on the 84th percentile load, obtained using
bootstrap techniques, are summarized in Table 5.13. The relative confidence bounds
(averaged over 10 sets) from 50 simulations, are less than 10% for both environmental
states, with linear and nonlinear waves. Earlier in Chapter 3, we suggested that a
given number of simulations is adequate to have a stable short-term distribution if
such a relative confidence bound is less than 15%, which is the case here. Moreover,
we see in Table 5.13 that the relative confidence bound decreases significantly when
all 500 simulations are used in a similar bootstrapping exercise. Besides, the 90% rel-
ative confidence bounds presented in Table 5.13, which are estimated using bootstrap
resamplings, are all smaller than the 82% confidence bounds shown in Table 5.12,
185
Table 5.14: Extrapolated loads (10-min maximum fore-aft tower bending moment at
mudline) from 50 simulations.
which were estimated from non-repeating sets of real data (albeit only ten sets). This
is also expected as bootstrap resamplings (with replacement) give smaller confidence
bounds than real non-repeating random samples [19]. In summary, the stable esti-
mates of a rare fractile and associated acceptable confidence bounds seen here suggest
that 50 simulations are adequate for extrapolation. We now discuss extrapolated 20
year tower loads for the offshore wind turbine being studied.
186
predicted loads (=100×(Max - Min)/Median) is less than 15% for all cases. As was
the case for the range on 84th percentile load in Table5.12, the computed range in
Table 5.14 is larger for nonlinear waves, indicating greater variability with nonlinear
waves. Also note that the extrapolated loads are very close to the maximum observed
load (over 500 simulations) for V = 18 m/s, Hs = 7.5 m; this is because the largest
load extreme from 500 simulations is associated with an exceedance probability of
1/501 or 1.996 × 10−3 , which is close to the desired (1 - p3 ) exceedance probability of
2.03 × 10−3 .
Linear Nonlinear
Wave Elevation (m)
−5
0 100 200 300 400 500 600
Wind Speed (m/s)
30
20
10
0 100 200 300 400 500 600
Fore−aft Tower Base Momemt (MN−m)
100
50
0
0 100 200 300 400 500 600
time (sec)
Figure 5.13: Ten-minute time series of the wave elevation, wind speed, and fore-aft
tower bending moment showing almost concurrent occurrences of the maximum wave
height and the maximum tower load for V = 18 m/s, Hs = 7.5 m.
187
Linear waves Linear waves
120 120
44 well correlated events 456 less correlated events
Maximum Load (MN−m)
100 100
90 90
80 80
70 70
6 8 10 12 6 8 10 12
Maximum Wave Height (m) Maximum Wave Height (m)
(a) (b)
100 100
90 90
80 80
70 70
6 8 10 12 6 8 10 12
Maximum Wave Height (m) Maximum Wave Height (m)
(c) (d)
Figure 5.14: Scatter plots showing whether the occurrences of the maximum fore-aft
tower bending moment in ten minutes and the occurrence of the maximum instanta-
neous wave height are close (i.e., well correlated) or far apart (i.e., less correlated).
Plots are for V = 16 m/s and Hs = 5.5 m.
To further understand how waves influence loads for the two environmental
states studied, we investigate whether the occurrence of the maximum fore-aft tower
bending moment in a ten-minute time series and the occurrence of the maximum
instantaneous wave height in that same time series are close (i.e., well correlated)
or far apart (i.e., less correlated) in time. A close occurrence would indicate that
188
Linear waves Linear waves
100 100
80 80
5 10 15 20 5 10 15 20
Maximum Wave Height (m) Maximum Wave Height (m)
(a) (b)
120 120
100 100
80 80
5 10 15 20 5 10 15 20
Maximum Wave Height (m) Maximum Wave Height (m)
(c) (d)
Figure 5.15: Scatter plots showing whether the occurrences of the maximum fore-aft
tower bending moment in ten minutes and the occurrence of the maximum instanta-
neous wave height are close (i.e., well correlated) or far apart (i.e., less correlated).
Plots are for V = 18 m/s and Hs = 7.5 m.
large waves might be causing large loads, for example in Fig. 5.13. We choose the
time separation (to distinguish well correlated from less correlated) to be equal to the
peak spectral period of the waves. Scatter plots of well and less correlated events in
each set are shown in Figs. 5.14 and 5.15 for V = 16 m/s and Hs = 5.5 m and for
V = 18 m/s and Hs = 7.5 m, respectively; the total number of such events are also
189
indicated on each of the plots. It is seen that for V = 16 m/s and Hs = 5.5 m, the
fraction of well correlated events is about 9% and 14% for linear and nonlinear waves,
respectively. For V = 18 m/s and Hs = 7.5 m, the fraction of well correlated events
is higher — about 21% and 44% for linear and nonlinear waves, respectively. In both
cases, the fraction of well correlated events is seen to get larger for nonlinear waves.
For V = 18 m/s and Hs = 7.5 m, this fraction is almost doubled, which is consistent
with our earlier observation that nonlinear waves for this environmental state result
in larger wave heights and, consequently, larger loads. For V = 18 m/s and Hs = 7.5
m, it can also be seen from the scatter plots (Fig. 5.15) that larger loads occur with
nonlinear waves, compared to those with linear waves. Furthermore, the scatter plot
shows strong correlation between maximum wave heights and maximum loads for the
nonlinear wave model only for V = 18 m/ sand Hs = 7.5 m. This clearly indicates
that nonlinear waves, especially for large wave heights, have a clear influence on loads.
Note that this is not the case for linear waves.
Our objective in this study was to investigate long-term loads for a utility-
scale 5MW offshore wind turbine sited in 20 meters of water. Our focus was on the
190
fore-aft tower bending moment at the mudline. We presented the theory for modeling
nonlinear (second-order) irregular waves, which is more appropriate than the use of
linear waves for shallow water depths, and have incorporated this nonlinear wave
model in a time-domain simulator that performs aero-servo-hydro-elastic analysis of
wind turbines. We have studied the influence of modeling nonlinear waves on loads
for a monopile support structure (a cylinder of 6 m diameter) of the turbine and
compared the loads derived to those based on the more conventional linear wave
theory.
We first focussed on the environmental state (mean wind speed of 16 m/s and
significant wave height of 5.5 m) that governed the 20-year load for this turbine, when
linear irregular waves were used. For this governing environmental state, we compared
loads due to linear and nonlinear waves on the monopile, both by a static analysis on
a rigid stand-alone monopile and by a dynamic analysis when the monopile is part of
the turbine.
• Loads due to waves alone are larger when nonlinear waves are modeled (com-
pared to when linear waves are modeled); the difference in overall loads dis-
appears when wind is also included as wind load effects on the turbine rotor
dominate waves for this environmental state.
• The wave steepness, for a given significant wave height, influences the hydrody-
namic loads such that monopile loads are largest at a higher steepness, while the
191
difference in loads computed from linear and nonlinear wave models decreases
with increasing steepness.
In summary, findings from this study suggest that nonlinear irregular waves
can have an important influence on the hydrodynamic loads experienced by an offshore
wind turbine support structure. Therefore, it is important that nonlinear waves be
considered when predicting long-term loads in evaluating ultimate limit states for
offshore wind turbine design.
192
Chapter 6
6.2 Conclusions
193
6.2.1 Offshore Environment at Blyth and its Effect on Offshore Wind
Turbine Support Structure Loads
Field data, when available, can be valuable in order to understand the envi-
ronment and its effects on the turbine response at an offshore site. Our objective in
Chapter 2 was to derive long-term loads for a 2MW offshore wind turbine sited in 9
meters of water at the Blyth site in the United Kingdom. Field data were available
on wind speed, wave height, and mudline bending moment in two formats, namely as
summary data (ten-minute statistics) and as campaign data (ten-minute time series).
The data were broadly classified into winds from the sea and winds from the shore.
Furthermore, the data were divided into 30◦ wind-direction sectors. Short-term load
distributions given environmental random variables were established using ten-minute
global maxima. Integration of these distributions with the relative likelihood of differ-
ent wind speed and wave height combinations allowed derivation of turbine long-term
loads.
• Winds from the shore govern the long-term loads for the instrumented wind
turbine at Blyth. A detailed analysis of the available time series showed that
winds from the shore result in larger loads, compared to winds from the sea, due
to: (1) the larger turbulence intensity associated with winds from the shore; and
(2) the relatively greater amount of energy in the waves at the tower natural
frequency (verified in load time series associated with winds from the shore
compared to those associated with winds from the sea).
194
• Uncertainty bounds on predicted long-term loads, derived using non-parametric
bootstrap techniques, can be quite large, and these were generally larger for
winds from the shore. These large uncertainty bounds result from the relatively
large variability in estimates of the parameters (as obtained from bootstrap
resamplings) used for the short-term load distributions. Furthermore, the width
of these uncertainty bounds was sensitive to the method used to estimate the
parameters of the short-term load distribution.
Using the LE3 data for loads on a model of a 5MW onshore wind turbine, we
have addressed several questions related to the theory and practical implementation
of statistical loads extrapolation. In particular, we focussed on establishing robust
short-term distributions, as they determine the accuracy of long-term loads predicted
using extrapolation, for a given site. This was our focus in Chapter 3. First, we
emphasized the need to understand which wind speeds tend to cause the largest
loads of different types. This information is useful for determining where greater
simulation effort is needed.
We compared the use of global and block maxima, two of several methods
that may be used to extract extremes from a time series, and their differences in
the context of statistical loads extrapolation. Issues related to the independence
of block maxima were addressed using a statistical test of independence. Ignoring
unimportant wind speed bins, it was found that block sizes of around 10-15 seconds
for four different load types (out-of-plane and in-plane bending moment at a blade
root, fore-aft tower bending moment at base, and out-of-plane blade tip deflection)
yielded independent block maxima when checked at the 1% significant level based
195
on the statistical test. Finally, though, with the LE3 data, it was demonstrated that
there is no advantage gained from using block maxima over global maxima when
short-term loads are estimated. Moreover, even when block sizes are very small so
as to clearly exhibit dependence, ignoring this leads to small error in estimation of
short-term loads.
196
wind as was the case for the onshore turbine studied; the procedure to establish such
short-term distributions for offshore turbines, however, remains virtually the same.
197
• Long-term load predictions based on the global maxima, block maxima, and
POT methods were found to be close as long as distribution tails were reliable
and well defined.
• The design wind speed governing the long-term loads was not the rated wind
speed (as is often the case) but was a wind speed above rated once load vari-
ability was taken into consideration.
• A chief source of load variability resulted from blade pitch control actions that
caused large loads such that the tails of the short-term load distribution were
not reliable unless a large number of simulations was performed.
These conclusions, while particular to the turbine model studied, are useful
to consider for any simulation-based exercise that seeks to predict long-term turbine
loads. This study also suggests that the effect of control actions on extreme loads
needs careful study; in particular it is of interest to investigate methods to account
for variability in loads due to control actions since they can alter the tails of the load
distributions in ways different than would be the case for loads not affected by control
actions.
198
the influence of modeling nonlinear waves on loads (focussing on the fore-aft tower
bending moment at the mudline) for a monopile support structure (a cylinder of 6
m diameter) of the turbine and compared the load statistics to those based on more
conventional linear wave theory.
• Loads due to waves alone were larger when nonlinear waves were modeled (com-
pared to when linear waves were modeled); this difference disappeared when
wind was also included since wind load effects on the turbine rotor dominated
waves for this environmental state.
• This monopile was dominated by inertia forces; the influence of nonlinear waves
would be even more pronounced on structures dominated by drag forces.
• The wave steepness, for a given significant wave height, influenced hydrody-
namic loads such that monopile loads were largest at a higher steepness, while
the difference in loads computed from linear and nonlinear wave models de-
creased with increasing steepness.
199
of 7.5 m). It was found that loads due to nonlinear waves, even when taking into
account the effect of wind loads on the turbine rotor, could be significantly larger than
those due to linear waves especially for drag-dominated support structures. Based
on simulations, we estimated empirical short-term distributions of the fore-aft tower
bending moment at the mudline due to linear and nonlinear irregular waves; it was
seen that, with nonlinear waves, loads corresponding to rare fractile levels, which are
of interest in predicting long-term loads using extrapolation, were significantly larger
than those predicted with linear waves.
In summary, findings from this study suggest that nonlinear irregular waves
can have an important influence on the hydrodynamic loads experienced by an offshore
wind turbine support structure. Therefore, it is important that nonlinear waves are
considered when predicting long-term loads in evaluating ultimate limit states for
offshore wind turbine design.
200
linear irregular waves on such long-term loads for turbine support structures. Based
on insights gained from this research, we recommend the following areas for further
research.
• The use of efficient extrapolation methods, such as the inverse FORM, in pre-
diction of long-term loads for offshore and onshore wind turbines for states other
than operating—namely, start-up, shut-down, and parked conditions—should
be investigated. Furthermore, the prediction of long-term fatigue loads, which
can often drive the overall design, using extrapolation procedures discussed in
this dissertation should be studied.
• The influence of nonlinear irregular wave models on fatigue loads for a turbine
support structure should be studied. In addition, support structures other than
a monopile, such as tripod and jacket structures, should be considered while
comparing different wave models.
• The effect of modeling the flexibility of foundations in ultimate and fatigue limit
states for wind turbine support structures should to be studied.
201
Bibliography
[1] American Wind Energy Association. AWEA 2007 Market Report. Washington,
DC, 2008.
[4] G. Birkhoff and J. Kotik. Fourier Analysis of Wave Trains. Proceedings, Sym-
posium on Gravity Waves, National Bureau of Standards, Circular No. 521,
Washington, D.C., 1951.
[7] E. A. Bossanyi. GH Bladed Version 3.6 User Manual. Technical Report Doc-
ument 282/BR/010 Issue 12, Garrad Hassan and Partners Limited, England,
2003.
202
[8] B. B. Brabson and J. P. Palutikof. Tests of the Generalized Pareto Distri-
bution for Predicting Extreme Wind Speeds. Journal of Applied Meteorology,
39(9):1627–1640, 2000.
[10] T. R. Camp et al. Design Methods for Offshore Wind Turbines at Exposed
Sites. Final Report of the OWTES Project, Garrad Hassan and Partners Ltd.,
Bristol, UK, 2003.
[16] Det Norske Veritas. Design of Offshore Wind Turbine Structures, Offshore
Standard, DNV-OS-J101, 2007.
203
[17] Det Norske Veritas. Environmental Conditions and Environmental Loads, Rec-
ommended Practice, DNV-RP-C205, 2007.
[19] B. Efron and R. J. Tibshirani. Statistical Data Analysis in the Computer Age.
Science, 253(5018):390–395, 1991.
[22] European Wind Energy Association. Delivering Offshore Wind Power in Eu-
rope — Policy Recommendations for Large-scale Deployment of Offshore Wind
Power in Europe by 2020. Brussels, Belgium, 2007.
204
[25] G. Z. Forristall. Irregular Wave Kinematics from a Kinematic Boundary Con-
dition Fit (KBCF). Applied Ocean Research, 7(4):202–212, 1985.
[28] M. Gilli and E. Kellezi. An Application of Extreme Value Theory for Measuring
Financial Risk. Computational Economics, 27:207–228, 2006.
[33] A. R. Henderson (ed.). Design Methods for Offshore Wind Turbines at Ex-
posed Sites - Hydrodynamic Loading on Offshore Wind Turbines. Technical
report, Delft University of Technology, Section Wind Energy, Stevinweg, The
Netherlands, 2003.
205
[34] W. Hoeffding. A Non-parametric Test of Independence. The Annals of Math-
ematical Statistics, 19(4):546–557, 1948.
[38] S.-L. J. Hu. Nonlinear Random Water Wave. In Computational Stochastic Me-
chanics, A.H.-D. Chang and C.Y.Yang, eds., Computational Mechanics Publi-
cation, Southampton, UK, 519-544, 1993.
[40] R. T. Hudspeth. Prediction of Wave Forces from Nonlinear Random Sea Sim-
ulations. Ph.D. Dissertation, University of Florida, 1974.
[41] R.T. Hudspeth and M.-C. Chen. Digital Simulation of Nonlinear Random
Waves. Journal of Waterway, Port, Coastal and Ocean Division, ASCE,
105(WW1):67–85, 1979.
206
[43] International Electrotechnical Commission. Wind Turbines — Part 3: Design
Requirements for Offshore Wind Turbines, IEC-61400-3, Edition 1.0, TC88
WG3 Committee Draft, 2007.
[44] A. K. Jha. Nonlinear Stochastic Models for Ocean Wave Loads and Responses
of Offshore Structures and Vessels. Ph.D. Dissertation, Stanford University,
1997.
[45] B. J. Jonkman and M. L. Buhl Jr. Turbsim User’s Guide, Revised February
1, 2007 for Version 1.21. Technical Report NREL/TP-500-41136, National
Renewable Energy Laboratory, Golden, CO, 2007.
[47] J. M. Jonkman and M. L. Buhl Jr. FAST User’s Guide. Technical Report
NREL/EL-500-38230, National Renewable Energy Laboratory, Golden, CO,
2005.
207
[50] N. D. Kelley. Full Vector (3-D) Inflow Simulation in Natural and Wind Farm
Environments using an Expanded Version of the SNLWIND (Veers) Turbulence
Code. Technical Report NREL/TP-442-5225, National Renewable Energy Lab-
oratory, Golden, CO, 1992.
[52] D. J. Laino and A. C. Hansen. User’s Guide to the Wind Turbine Aerodynamics
Computer Software AeoDyn. Technical report, Windward Engineering, LC, Salt
Lake City, Utah, 2002.
[54] T. J. Larsen and A. M. Hansen. How 2 HAWC2, The User’s Manual. Technical
Report Risø-R-1597, ver 3.1, Risø National Laboratory, Technical University of
Denmark, Rosklide, Denmark, 2007.
208
[58] M.S. Longuet-Higgins. The Effect of Non-linearities on Statistical Distributions
in the Theory of Sea Waves. Journal of Fluid Mechanics, 17(3):459–480, 1963.
[61] P. Lynett and P. L.-F. Liu. A Two-layer Approach to Wave Modelling. Pro-
ceedings of the Royal Society of London, Series A, 460:2637–2669, 2004.
[63] P. A. Madsen, H. B. Bingham, and H. Liu. A new Boussinesq method for Fully
Nonlinear Waves from Shallow to Deep Water. Journal of Fluid Mechanics,
462:1–30, 2002.
[64] P. H. Madsen, K. Pierce, and M. Buhl. Predicting Ultimate Loads for Wind
Turbine Design. In ASME Wind Energy Symposium, AIAA-99-0069, Reno, NV,
1999.
[65] R. E. Melchers. Structural Reliability Analysis and Prediction. John Wiley, New
York, 1999.
[66] Minerals Management Service. Technology White Paper on Wind Energy Po-
tential on the U.S. Outer Continental Shelf. Renewable Energy and Alternate
Use Program, U.S. Department of the Interior, Washington, DC, 2006.
209
[67] P. J. Moriarty. Database for Validation of Design Load Extrapolation Tech-
niques. Wind Energy, to be published, 2008.
[70] W. Musial and S. Butterfield. Future for Offshore Wind Energy in the United
States. Technical Report NREL/CP-500-36313, National Renewable Energy
Laboratory, Golden, CO, 2004.
[73] L. D. Nelson. Analysis of Wind Turbine Loads using Inflow and Structural
Response Data from Field Measurements. M.S. Thesis, The University of Texas
at Austin, 2003.
[74] E. Norton. Investigation into IEC Offshore Draft Standard Design Load Case
1.1. RECOFF Document 59, Recommendations for Design of Offshore Wind
Turbines, Denmark, 2004.
210
[75] O. Nwogu. Alternative form of Boussinesq Equations for Nearshore Wave
Propagation. Journal of Waterway, Port, Coastal and Ocean Engineering,
119(6):618–638, 1993.
[77] S. Øye. FLEX5 Users Manual. Technical report, Technical University of Den-
mark, Lyngby, Denmark, 1999.
[79] W. J. Pierson, Jr. and P. Holmes. Irregular Wave Forces on a Pile. Journal of
Waterways and Harbors, 91(WW4):1–10, 1965.
[80] E. Powers. Class Notes for EE 381L: Higher-order Statistical Signal Processing.
The University of Texas at Austin, 2006.
211
[84] A. A. Ramos. Extreme Value Theory and the Solar Cycle. Astronomy and
Astrophysics, 472(1):293–298, 2007.
[89] K. Saranyasoontorn and L. Manuel. Design Loads for Wind Turbines using
the Environmental Contour Method. Journal of Solar Energy Engineering,
Transactions of the ASME, 128(4):554–561, 2006.
212
Wave Forces. Technical Report Ocean Engineering Report No. 20, University
of Delaware, Newark, DE, 1979.
[94] M. Shinozuka and C. M. Jan. Digital Simulation of Random Processes and its
Applications. Journal of Sound and Vibration, 25(1):111–128, 1972.
[97] B. Sweetman and S.R. Winterstein. Second-order random Ocean Waves: Pre-
diction of Temporal and Spatial Variation from Fixed and Moving References,
The Routine WaveMaker. Technical Report RMS-37, Civil Engineering Depart-
ment, Stanford University, May 1999.
213
Conference on Offshore Mechanics and Arctic Engineering, OMAE1995, pages
313–319, Volume 1A, pp. 313-319, Copenhagen, Denmark, 1995.
[101] J. van der Tempel and D.-P. Molenaar. Wind Turbine Structural Dynamics - A
Review of the Principles for Modern Power Generation, Onshore and Offshore.
Wind Engineering, 26(4):211–220, 2002.
[102] The Specialist Committee on Environmental Modeling. Final Report and Rec-
ommendations to the 22nd ITTC. Technical report, International Towing Tank
Conference, Seoul, 1999.
214
[108] J. Zhang, L. Chen, M. Ye, and R. E. Randall. Hybrid Wave Model for Uni-
directional Irregular Waves — Part I. Theory and Numerical Scheme. Applied
Ocean Research, 18:77–92, 1996.
215
Vita
Puneet Agarwal was born in Lucknow, India on 25 April 1977, the son of Mr.
Chandra Kishore Agarwal and Mrs. Savitri Devi Agarwal. He received the degree
of Bachelor of Technology in Civil Engineering from Indian Institute of Technology
Kanpur in 2000, and the degree of Master of Science in Civil Engineering from Uni-
versity of Connecticut in 2004. He entered the doctoral program in the Department
of Civil, Architectural and Environmental Engineering at the University of Texas at
Austin in 2004.
† A
LT EX is a document preparation system developed by Leslie Lamport as a special version of
Donald Knuth’s TEX Program.
216