Seismic Load Evaluation of Wind Turbine Support Structures with

Consideration of Uncertainty in Response Spectrum and Higher Modes

T. Ishiharaa, G. Takamotob, M. W. Sarwarc
a
ishihara@bridge.t.u-tokyo.ac.jp b takamoto@bridge.t.u-tokyo.ac.jp c sarwar@bridge.t.u-tokyo.ac.jp
a,b,c
Department of Civil Engineering, The University of Tokyo.
7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656 JAPAN

Summary
Rapid development of wind energy in seismically active regions like Japan requires evaluations of
design seismic load for support structures to ensure structural integrity. This paper presents a
rational approach to determine design response spectrum for structures with low damping ratio
such as wind turbine towers. A modified damping correction factor for the design response
spectrum of wind turbine towers is proposed to incorporate uncertainty in seismic loads. Proposed
correction coefficient of the design response spectrum is defined as a function of percentile quantile
of the seismic response distribution, which is determined by code calibration method against the
current seismic design loads. This study also presents a simplified formula for the seismic load
evaluation of wind turbine towers based on response spectrum analysis. First three tower modes
are found to have significant contributions to the seismic load of wind turbine towers and expected
maximum seismic load is obtained by combining three modal responses using SRSS method.
Expressions for estimation of first three mode shapes and natural periods are proposed to facilitate
the routine design implementation. Finally the accuracy of proposed design response spectrum and
seismic load formula is evaluated through comparison of results with time history analysis for wind
turbines of capacities from 400kW to 2MW.

1. Introduction

Estimation of seismic response of wind turbines becomes of great importance when wind farms are
designed and developed in seismically active regions. Being a seismically active region, Japan has
strict regulations for designing and assessing the safety of wind farms. This require use of design
formula based on response spectrum method and time domain analysis for low and high structures
respectively. Current practice of estimating design loads based on the response spectrum method
[1]
encounter two problems when applied to the wind turbine support structures . These support
structures are extremely low damped and experience a wide range of frequencies when subjected
to seismic activities. Response spectrum for such low damped structures show excessive
fluctuation and such uncertainty in response spectrum can not be captured by existing models of
[2] [3]
the damping correction factors defined in Eurocode and BSL . In addition, use of the simplified
[4]
SDOF model suggested by IEC results in linear vertical load profiles. However, vertical distribution
[1]
of the seismic loads is found to be largely affected by the higher modes of wind turbines.
Therefore simplified but accurate analysis method to estimate design load profiles is desired.
In this research, a modified damping correction factor that accounts for uncertainty in
response spectrum, and modal participation functions encompassing complex vertical distribution of
seismic loads are proposed. The accuracy and reliability of the proposed method for evaluation of
seismic design loads is examined against time history analysis and current design codes.

2. Seismic Load Estimation

Generally the seismic design loads are estimated by two methods, time history analysis and the
response spectrum method. In time history analysis, equation of motion is solved for different
earthquake waves to obtain loads, shear and bending moments, acting on the wind turbine support

a
Professor bGraduate Student cResearcher
Phone: +81 3 5841 1145 Fax: +81 3 5841 1147
-1-
structure. This method has inherent advantage of including the effect of higher modes and
geometrical non-linearity on the structural response when structure is subjected to dynamic
excitation. However, this method is relatively time consuming and requires certain number of
seismic waves to account for variability of the analysis results. However, response spectrum
method requires only natural period, mode shape and mass distribution of the structure to calculate
maximum seismic loads.
This section describes basic assumptions used in the seismic load analysis, introduces
equation of motion for time history analysis, response spectrum analysis method along with
acceleration spectrum and, finally details of the wind turbine models used in this study are
discussed.

2.1 Analysis assumptions

Basic assumptions for estimation of seismic loads acting on the wind turbine support structures are
listed below:
a. Fixed support model is used for eigen value and dynamic analysis of the wind turbine
support structure. Seismic waves at the tower base level are used.
b. The wind turbine tower is modeled as multi-degree of freedom system. Rotor and
nacelle masses are modeled as a lumped mass situated at hub height.
c. To determine seismic loads under operating conditions, seismic load under parked
[5]
conditions is combined with the operational wind loads
d. Additional loads caused by geometrical non-linearity such as P-∆ effect and twisting at
[7]
top of the tower are considered separately.

2.2 Time History Analysis

In time history analysis, seismic loads acting on the wind turbine are estimated by introducing the
synthetic earthquake waves at the tower base. These waves are generated to closely match the
target spectrum defined by equation (4). In this study, wind turbine tower is modeled as a multi-
degree of freedom (MDOF) system, as shown in Figure 1, where mass of tower members are
lumped at nodes and beam element model is used to model stiffness and damping of the respective
members. Under parked conditions, the equation of motion for MDOF system subjected to ground
acceleration
x
g is:

[m]{
x
} + [c]{x
} + [k ]{x} = −[m]{e}
x
g (1)

where {
x
}, { x
}, { x} are the acceleration, velocity and displacement vectors respectively.
[m ] , [c ] and [k ] are the mass, damping and stiffness matrices. {e} is a unit vector.

2.3 Response Spectrum Method

Response spectrum method is based on acceleration response spectrum of SDOF system that is
used to determine maximum response of MDOF system. Following is the eqation of motion for jth
mode of a MDOF system,

q

j + 2ζ jω j q j + ω 2j q j = −γ j
x
g (2)

， ，
where ω j is the natural angular frequency ζ j is the damping ratio γ j is the modal participation
factor. Maximum force for each mode depends upon the modal participation factor, mode shapes
and acceleration of response spectrum corresponding to natural period of respective modes. For
example, the maximum shear force corresponding to j-th mode of the i-th node can be estimated as
follows:

Fij = γ j X ij S a (T j , ζ )mi (3)

-2-
where X ij is the j-th mode shape, S a (T , ζ ) is the amplitude of acceleration response spectrum
corresponding to the natural period T and damping ratio ζ . To estimate maximum design loads,
first load for each mode of the MDOF system is determined using equation (3) and then loads for all
modes are combined. As response spectrum method uses dynamic characteristics such as mode
shapes and natural period of the structure, a prior knowledge of these characteristics is required to
estimate the seismic loads. Therefore, to facilitate the load estimation procedure, this study
proposes a model to estimate these characteristics without eigen value analysis.

Seismic input

SDOF MDOF

Fig. 1 Fixed support models used for load analysis

2.4 Acceleration Response Spectrum

[2]
The design acceleration response spectrum (Sa) used to determine seismic loads is defined as ,

 T
a 0 ⋅ S ⋅ {1 + ⋅ ( β 0 ⋅ Fζ − 1)} (0 ≤ T ≤ TB )
 T B

S a (T , ζ ) = a 0 ⋅ S ⋅ Fζ ⋅ β 0 (TB ≤ T ≤ TC ) (4)

a 0 ⋅ S ⋅ Fζ ⋅ β 0  TC 
⋅  (TC ≤ T )
 T 
α
 7 
Fζ (ζ ) =  
 (5)
 2 + 100ζ 

Where a0 is the peak ground acceleration at the engineering bed rock for a given return period, S is
the soil amplification factor, Fζ is the damping correction factor and β0 is acceleration response
magnification ratio for the region where acceleration response becomes constant, TB and Tc
defines the range of constant spectral acceleration. Parameters to define a design acceleration
response spectrum with a return period of 500 years[5] for hard ground strata, which is defined as
Soil Type 1, are listed in Table 1.

Table 1 Response spectrum parameters used in this study

a0 (m/s 2 ) S β0 TB (s) TC (s)
3.2 1.5 2.5 0.16 0.576

2.5 Description of wind turbine models

Previous study has shown that the ratio of mass of nacelle and blades to the total mass of the wind
[5]
turbine system remains almost constant regardless of the rated power of wind turbines as shown

-3-
in Table 2. Therefore, in this study, two wind turbines with rated power of 400kW and 2MW are
selected as typical examples for investigations. Details of basic parameters of these wind turbines
are summarized in Table 3.

Table 2 Mass ratio of different wind turbines Table 3 Details of the wind turbine models
Rated Power Mass Description Unit Model 1 Model 2
Model
(kW) Ratio Rated Power kW 400 2000
1 100 0.44 Rotor diameter m 31 80
2 400 0.46 Rotor tilt deg 5 5
3 500 0.45 Tower height m 35 67
4 1000 0.41 Hub height m 36 67
5 1500 0.45 Blade mass kg 1100 6800
6 2000 0.44

3. Proposal for estimation of seismic loads

Application of the response spectrum method to estimate seismic loads of the wind turbine support
structure is hindered by two problems. Wind turbine support structures are significantly low damped,
and response spectrum for such low damped structures show excessive fluctuation. Such
uncertainty in response of these structures can not be captured by existing models of the damping
[2]
correction factors defined in Eurocode . Therefore, formulation of new damping correction factor is
needed to incorporate extreme fluctuation in spectral acceleration to establish a reliable design
spectrum for wind turbines. In addition, vertical distribution of the seismic loads is found to be
[1]
largely affected by the higher modes of wind turbines. However, previous estimation of vertical
[4]
load profile is based on the simplified SDOF model suggested by IEC that fails to capture the non-
linearity of vertical load profiles. Therefore it is necessary to include higher modes in the load
estimation procedure.
In this section, a damping correction factor that account for uncertainty in the response
spectrum and modal participation function for higher modes are proposed for accurate estimation of
the design loads.

3.1 Model for damping correction factor

To account for excessive fluctuations in the response spectrum of low damped systems, damping
correction factor is proposed as a function of spectral uncertainty, natural period and damping ratio
so that,

α
 7 
Fζ (ζ , T , γ ) =   , α = f (T , γ ) (6)
 2 + 100ζ 

where T is natural period function, ζ is damping ratio and γ is the quantile value for desired reliability
level. Previously Eurocode defined damping correction factor as a function of damping ratio only
[2]
and a constant value of 0.5 was used for the exponent α . However, when exponent α is defined
as a function of time period, f(T), it corresponds to a damping correction factor that considers the
[7]
natural period of the structure . In addition to damping ratio and natural period of the structure,
proposed damping correction factor also includes uncertainty of the response spectrum.
To establish the proposed damping correction factor, first a set of seismic waves is
generated to statistically evaluate the uncertainty in response spectrum of SDOF system. Then
relation for exponent α was established by data fitting to different quantiles of the response
spectrum. Following describes the procedure in detail.

-4-
a. Probability distribution of uncertainty in Response Spectrum
A set of 35 seismic waves, 5 with observed phase and 30 with random phase, was used to evaluate
excessive fluctuation in the acceleration response spectrum for range of damping ratios ,i.e., from
0.5% to 5%. Figure 2 shows acceleration response spectra for damping ratios of 0.5% and 5% that
correspond to wind turbine structures and buildings respectively. It can be observed that at low
damping ratio of 0.5%, in case of wind turbine support structures, excessive fluctuations in the
spectral acceleration occur.

1
5000

Hachinohe

Cumulative Relative Frequency

]l4000 Taft
0.8

ag
[ Elcentro
no3000 Kobe 0.6
it
ar
el2000 ζ=0.5%
ec 0.4
cA
1000 Section I A
IB 0.2 Section I B
IA Section I C
0
IC
ζ=5% Log Normal Dist
0.1 1
0
Time Period T[s] 1000 10 4
Fig. 2 Sections of acceleration response spectrums Fig. 3 CRF of spectral acceleration for
for statistical investigation each section

To determine probability distribution that represents uncertainty involved, response spectrum is

divided into three sections so that 0.05<T<TB refers to section IA，TB≤T≤TC refers to section IB
and TC<T<5 refers to section IC. Sections IA and Ic that define the non-linear regions of the
response spectrum are divided into ten sub-sections I A(i ) ・ I C(i ) (i = 1 ~ 10) , whereas section IB was
considered as a single section to calculate the statistical properties such as mean and standard
deviation of acceleration response. Figure 3 shows a cumulative relative frequency in the sections
I A( 5) ・ I B ・ I C( 5) . Also cumulative frequencies of lognormal distribution function derived from mean and
standard deviation for each interval are drawn as solid lines. Log normal distribution is found to
have well defined the uncertainty in all sections of the acceleration response spectrum as shown in
above Fig. 3.

b. Identification of exponent α in damping correction factor

It is now possible to define the percent quantile γ of acceleration for a desired reliability level by
modeling the uncertainty of acceleration response with logarithm normal distribution function. Three
quantile values of 20%, 50% and 80% are used to investigate exponent α of the damping correction
factor. Percentile quantile value of acceleration spectrum is used for Sa(T,ζ) of equation (4) to
calculate damping correction factor Fζ and then exponent α is calculated using equation (5). The
exponent α is found to have a linear relation with percent quantile γ and natural period T, as shown
in Fig. 4 and Fig. 5, that leads to following:

α = f (T , γ ) = −0.07T + 0.7γ + 0.5 (7)

The acceleration response spectrum calculated by proposed formula agrees well with the calculated
results for all quantile values γ of response acceleration well. Introduction of natural period T of

-5-
structure has lead to accurate estimation of response spectrum in the long period regions. Also
uncertainty of the response spectrum can be incorporated by changing quantile value γ. However,
for response spectrum based on the damping correction factor defined in Eurocode, it is found that
it corresponds to 20% quantile values.

1.5
1.5
ζ=0.5% γ=0.2
γ=0.5
γ=0.8
1
1

α α
0.5
0.5
T=0.304
T=1.52
T=4.49
0
0 0 1 2 3 4 5
0 0.2 0.4 0.6 0.8 1
Quantile　γ Time Period　T[s]
Fig. 4 Variation of α with quantile γ Fig. 5 Variation of T with quantile γ

4000
20% Quantile(γ=0.2)
50% Quantile(γ=0.5)
]
l3000 80% Quantile(γ=0.8)
a
g Proposed Eq
[ Eurocode
n
o
i2000
t
a γ=0.8
r
e γ=0.5
l
e
c γ=0.2
c 1000
A

0
0.1 1
Time Period　T[s]
Fig. 6 Variation of response spectrum with γ-values

3.2 Vertical load distribution function considering higher modes

Seismic loads acting on wind turbine support structure can be calculated using load formula for
MDOF system presented in equations (8) and (9). But estimation of design loads for MDOF system,
shown in Fig 1, require prior knowledge of natural periods and modal participation function γjXij for
predominant modes. These structural characteristics are generally calculated by performing eigen
value analysis. However, this study proposes use of natural period ratio and polynomial expression
to estimate the natural period and modal participation function of higher modes. Then SRSS
method is used for super positioning of these modes to obtain the maximum design load for wind
turbines as shown in equation 10.

n n
Qij = ∑ Fkj = ∑ γ j X kj S a (T j , ζ )m k (8)
k =i k =i
n n
M ij = ∑ Fkj ( z k − z i ) = ∑ γ j X kj S a (T j , ζ ) m k ( z k − z i ) (9)
k =i k =i

n n
Qi = ∑Q j =1
ij
2
Mi = ∑M
j =1
ij
2
(10)

-6-
a. Modeling of modal participation functions and time period ratio
Mode shape of wind turbines are generally calculated using eigen value analysis. Previous research
[8] presented a polynomial expression as a function of height ratio for the modal participation
function of first mode γ1Xi1 of wind turbine support structure. This study proposes polynomial
equation as a function of height ratio to estimate the modal participation function for higher modes
of the wind turbine support structures as shown below:
k +1
 zi 
γ j X ij = ∑c
k =1
jk  
H 
(11)
where zi, H and cjk are the height of i-th node, hub height of wind turbine tower and coefficient of
polynomial respectively. The time period of higher modes are presented as a ratio of the natural
period of respective mode to that of the first mode, e.g., Tj /T1 is defined for period of the j-th mode.
Table 4 shows the coefficient of polynomial cjk and time period ratio Tj/T1 for the first three modes.
1

Table 3 Coefficients of polynomial and 0.8

time period ratio

Height Ratio
c jk 0.6

k
T j / T1
j 1 2 3
0.4
1 1.1 0 0 1
2 5.87 -6.00 0 0.127 1st Mode
0.2
2nd Mode
3 14.26 -38.20 24.00 0.043 3rd Mode
Equation
0
-1.2 -0.8 -0.4 0 0.4 0.8 1.2

Model participation Function( γ j Xij )

Fig. 7 Comparison between the proposed
& calculated γjXij

Fig. 7 shows comparison between eigen and proposed modal participation function for first three
modes of 400kW and 2MW wind turbine. From this figure, it is clear that modal participation function
of different sized wind turbines can be calculated accurately by the proposed polynomial ．
b. Estimation of seismic load for each mode and their combination
For load estimation by spectrum method, IEC61400 requires use of structural modes that account
[4]
for the required total modal mass of 85% . In this study, first three modes are found to satisfy this
criterion. First shear load profile for selected modes are calculated using the proposed equations of
modal participation function and time period ratio as shown in Fig 8. In case of 400kW turbine,
st
contribution of higher modes is smaller compared to 1 mode. However, a significant contribution by
nd rd
2 and 3 modes can be observed at the base of 2MW turbine.
1 1

400kW 2MW
0.8 0.8

o 0.6 o 0.6
tia i
t
R a
R
ht t
h
ig g
i
He 0.4 e 0.4
H

1st Mode 1st Mode

0.2 0.2 2nd Mode
2nd Mode
3rd Mode
3rd Mode

0 0
-400 -200 0 200 400 -1500 -1000 -500 0 500 1000 1500
Shear Force Qj[kN] Shear Force Qj[kN]

Fig. 8 Shear force profile of first three modes (γ = 0.5)

-7-
Since consecutive modes of wind turbine are well separated, SRSS method is used for
superposition of these modes as shown in equation(10). Figure 9 shows profile of shear force ratio
of current study by SRSS method and time history analysis along with previous work[8] based on
the 1st mode of vibration. It is clear that load profiles obtained from proposed equations could
capture the non-linearity of vertical profile and show good agreement with time history analysis
results.
1 1

0.8 0.8
Height Ratio

Height Ratio
0.6 0.6

0.4 0.4

0.2 0.2
Previous Model Previous Model
Proposed Model Proposed Model
Time History Analysis Time History Analysis
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Shear Force Ratio Moment Ratio

Fig.8 Vertical profiles of shear force and bending moment ratios

3.3 Determination of quantile γ for current reliability levels

In section 3.1, quantile γ was introduced to the damping correction factor that accounted for large
uncertainty of response spectrum. Therefore, it is necessary to determine suitable quantile value γ
[9]
for defining reliable design spectrum. Based on reliability theory, code calibration method is an
effective approach to ensure essentially similar reliability level as that of current design codes. In
[4]
Japan, BSL requires time history analysis of structures for atleast three earthquake waves to
obtain structural design certification. These waves include two local waves, such as Kobe and
Hachinohe that are onshore & offshore waves respectively, and one famous earthquakes like
Elcentro or Taft. In this study, four of the above mentioned seismic waves are used to determine
suitable quantile γ for determining the design response spectrum. The vertical load profiles obtained
by time history analysis of selected wind turbines are shown in Fig. 9. A γ-value of 0.7, i.e., 70%
quantile, is identified against the time history analysis results to obtain reliability level similar to that
established by the current design code.

400kW 2MW
1 1
Hachinohe Hachinohe
Taft Taft
0.8 Elcentro 0.8
Elcentro
Kobe Kobe
o γ =0.7 o γ=0.7
i
t 0.6 i
a t 0.6
R a
R
t
h t 2MW
g h 0.4
i 0.4 g
i
e e
H H
400kW
0.2 0.2

0 0
0 500 1000 15 00 2000 2500 0 2 10 4 4 10 4 6 10 4 8 10 4 1 10 5
Shear Force [kN] Shear Force [kN]

Fig. 9 Comparison of seismic load profiles: Time history analysis and proposed formula (γ = 0.5)

-8-
4. Conclusions

In this study, a modified damping correction factor is proposed that accounts for the excessive
fluctuation of response spectrum at low damping ratios and consider natural period of the wind
turbines. In addition, formula for analytical estimation of complex profile of seismic design loads are
presented that introduces contribution of higher modes, up to third mode, to the vertical load
distribution. Finally accuracy of proposed formula is verified against time history analysis and,
reliability level similar to that established by the current design code is demonstrated using code
calibration method.

References
1. Ishihara T, Sarwar MW. Numerical and Theoretical Study on Seismic Response of Wind
Turbines. Proc. of EWEC 2008: 2008.
2. Eurocode 8: Design of structure for earthquake resistance;Part1:Genersl rules, seismic actions
and rules for buildings,1998-1:2004.
3. BSL. The Building standard law of Japan; The building centre of Japan, 2004. (in Japanese)
4. IEC61400-1. Wind turbines. Part 1: Third edition.
5. JSCE. Guidelines for design of wind turbine support structures and foundations; Japanese
society of civil engineers, 2007. (in Japanese)
6. Nishimura I, Noda S, Takemura K, Ohno S, Tohdo M, WatanabeT. Response spectra for design
purpose of stiff structures on rock sites. Trans. SMiRT16 2001; P. No. 1133.
7. JSCE. Guidelines for design of wind turbine support structures and foundations; Japanese
society of civil engineers, 2010. (in Japanese)
8. Ishihara T, Zhu L, Binh LV. Earthquake load estimation method for the wind turbines under
operation and parked conditions. 29th Wind Engineering Symposium: Tokyo, 2007; 187-190
9. Hoshiya M, Ishii K. Reliability based design for structures; Kajima Press: Tokyo, 1986.

-9-